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- ABSTRACT 


The objective of this work is the investigation of vortical flows at high angles of at- 
tack using numerical techniques. , First step for a successful application of a numerical 
technique, such as finite difference or finite volume, is the generation of a computational 
mesh which can capture adequately and accurately the important physics of the flow. 
Therefore, the first part of this work deals with the grid generation over a double-delta 
Wing and the second part deals with tlie visualization of the computed flow field over the 





double-delta wing at different angles of attack. The surface geometry of the double-delta 
wing is defined algebraically. The developed surface grid generator provides flexibility 
in distributing the surface points ulong the axial and circumferential directions. The 
hyperbolic grid generation method is chosen for the field grid generation and both cv- 










lindrical and spherical grids are constructed. The computed low speed (M = 0.2) flow 

results at different angles of attack over the double-delta wing are visualized. Important 

flow characteristics of the leeward side flow field are discussed while the development 
: of vortex interaction. occurrence and progression of vortex breakdown as the angle of 

attack increases is demonstrated. The computed results at different fixed angles of at- 
" tack are presented. ae 
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THESIS DISCLAIMER 


The reader is cautioned that computer programs developed in this research may not 
have been exercised for all cases of interest. While every effort has been made, within 
the time available, to ensure that the programs are free of computational and logic er- 
rors, they cannot be considered validated. Any application of these programs without 
additional verification is at the risk of the user. 
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I. INTRODUCTION 


The objective of this work is the investigation of vortical flows over three- 
dimensional bodies at high incidence utilizing numerical methods. The advantage of 
numerical simulations compared with experiment is that they allow simultaneous obser- 
vation of all flow quantities of interest for the entire flow field. The disadvantage of 
numerical techniques is the accuracy limitations for simulations of flow fields over 
complex realistic configurations, even with the most efficient numerical schemes and fast 
computers. High Reynolds number turbulent flows of engineering interest can be fully 
simulated when all relevant scales are resolved. Resolution of all scales for complex high 
revnolds number flows over realistic configurations is beyond the capabilities of the 
present and next generation supercomputers. Common practice for the simulation of 
engineering flows is the use of various turbulence models to approximate the effect of 
the small scales which cannot be resolved. Error sources in numerical simulations are 
related to the discretization process, the order of accuracy of the numerical scheme and 
the turbulence modeling that is used. 

Nevertheless. Computational] Fluid Dynamics (CFD) allows investigation of various 
fluid flow phenomena that in the past was possible only in wind tunnels, water tunnels 
or actual flight testing. The advantage of being able to accurately capture the flow 
characteristics over complex configurations or even complete aircraft without endan- 
gcring life, 1.e. prelinunary flight testing. is readily apparent. Numerical solutions also 
enable to investigate and visualize the flow field characteristics from anv viewpoint or 
in as much detail as desired. With the ever increasing speed cost ratio of today’s com- 
puters. CFD techniques will be playing a more significant role facilitating aerodvnamic 
research and supplementing experimental investigations. Even though CFD and 
Navier-Stokes methods are not a new research tool, new and more efficient numerical 
techniques are evolving, while at the same time computers are becoming faster. Nu- 
merical prediction of steady flows over complete aircraft and comparison with flight data 
is already underway [Ref. 1]. In the near future CFD is expected to play a more active 
role in fluid dynamic research enabling simulation of complex unsteady flow regimes. 

In the past panel methods and vortex lattice methods were used in the analvsis of 
flows [Ref. 2]. These methods were insufficient for a detailed analysis of complex Nows 


such as vortical flows over bodies at incidence. The limitations of these methods are due 





to the potential flow assumption which is valid only for inviscid and irrotational flow. 
Viscous effects close to the surface for attached or mildly separated flow are obtained 
using Boundary Laver methods [Ref. 3]. The rotational compressible flow regime at high 
Reynolds numbers was investigated with the Euler «quations. Viscous effects become 
more important for flows at high angles of attack; therefore, the solution of the Navier- 
Stokes equations is required for this flow regime. 

In Chapter 2 the theoretical development of the compressible Navier-Stokes 
equations will be discussed. A finite difference algorithm used for the numerical solution 
of these equations will be presented. The numerical solution is performed on the finite 
number of points obtained after discretization of the flow domain. The procedure which 
vields this finite collection of points in the solution domain is known as grid generation. 
The quality of the solution depends directly on the smoothness of the grid and its ability 
to accurately represent flow gradients. Therefore, grid generation is an important part 
of the numerical solution. However, the numerical solution of the governing equations 
iS not the main objective of this research. The grid generation part. which is a necessarv 
stage before starting the numerical solution will be covered in full detail. Numerical 
solutions depend on the representation of the flow field by an orderly, finite collection 
of points. The process of obtaining three-dimensional grids involves first definition of 
an inner boundary. commonly known as the surface grid, before the subsequent gener- 
ation of the field grid can begin. 

The methods available for both the surface and field grid generation will be covered 
in Chapters 3 and 4, respectively. Developments in the area of grid generation have 
provided a kev to eliminate the problem of boundary shape definition [Ref. 4]. Finite 
difference grids can also be used to construct meshes that are suitable in finite element 
methods. The specific numerical method utilized in this research is the finite difference 
method. The finite difference method is one of the oldest numerical methods that can 
be utilized to obtain numerical solutions to differential equations. The application of 
this method is based on a Taylor series expansion and the definition of the derivative: 
most likely first developed by Euler in 1768 [Ref. 5: p. 167]. The algorithm used for the 
numerical integration utilizes a partially Mux-split numerical scheme with central differ- 
encing in the other two directions (Ref. 6]. 

The methods described above will be applied to a double-delta wing that has a strake 
with a sweep angle of 76° and a delta wing section with a sweep angle of 40°. Particular 
emphasis will be placed on the investigation of the vortical flow field at moderate to high 


angles of attack. Separated flow along the strake’s leading edge forms free shear lavers 

















which roll up to form vortex cores. This primary strake vortex generates an additional 
non-linear lift called vortex induced lift. A primary wing vortex also develops from the 
leading edge of the 40° swept delta wing. 

The mutual interaction of the strake and wing vortices and their interaction with the 
surface is an active area of current research. Investigation and prediction of the vortex 
breakdown that appears at higher angles of attack is also of high interest. The devel- 
opment of the leading edge vortex as well as breakdown are important phenomena that 
need to be fully understood. VWarious angles of attack, a = 10.0°, 19.0°, 22.4° are in- 
vestigated and compared with available experimental data. Understanding the leeward- 
side flow structure as well as breakdown are important phenomena that affect 
significantly today’s tactical and fighter aircraft effectiveness. 

Vortex breakdown is a transition of the vortex core from a jet-like flow to a wake- 
hke flow. Both swirl angle and adverse pressure gradient along the axial direction con- 
tribute to the breakdown of the vortex. Peckham and Atkinson first identified vortex 
breakdown when analvzing delta wings at high angles of attack [Ref. 7]. Research on 
vortex breakdown was continued by Elle, Lambourne and Brver, Harvey, Pritchard, 
Sarpkava, Hummel, Faler and Leibovitch, Payne and Nelson [Ref. 8,9,10,]1, 12 
.13,]4.15,16]. Studies then naturally progressed to more complicated bodies such as the 
double-delta wing where Brennenstuh] tested several wings 1n a Jow speed wind tunnel | 
and a water tunnel [Ref. 17]. The present study will attempt a comparison of the com- 
putational solution with the data obtained from wind tunnel testing done by 
Cunningham and Boer [Ref. 18]. This comparison along with discussions of the results 
that were developed during this research wiil be covered in Chapter 5. The closing 


chapter summarizes the conclusions and presents recommendations for further research. 














Il. THEORETICAL APPROACH 


The main objective of this work is the investigation of different techniques for the 
grid generation over complex three-dimensional bodies, and numerical flow visualization 
of the computed flowfields over bodies at high angles of attack. The flow field is ob- 
tained by the numerical solution of the Navier-Stokes equations. Fluid flow in the 
continuum flow regime includes most of the physical flows and is governed by the 
Navier-Stokes equations. The derivation of the Navier-Stokes equations is well known 
[Ref. 3: pp. 47-66]. Solutions of these equations are of interest in basic fluid mechanics 
research and for engineering applications. The solution of the Navier-Stokes equations 
is quite difficult due to their nonlinearity. Analytical closed form solutions of the 
Navier-Stokes equations can be obtained for only a few flow situations of simple ge- 
ometrical configurations and boundary conditions. Simplified forms of the Navier- 
Stokes equations, such as the boundary laver equations, can give satisfactory answers 
for many flows of practical interest. The inviscid form of the Navier-Stokes equations, 
commonly known as the Euler equations. can provide solutions for flows away from 
solid boundaries. However, complex flows such as vortical separated flows require the 
solution of the full Navier-Stokes equations. which can only be obtained by utilizing 
numerical techniques. The derivation of the Navier-Stokes equations is outlined in the 
following paragraphs. 


A. GOVERNING EQUATIONS 
1. The Continuity Equation 

For the derivation of the Navier-Stokes equations the fluid medium js consid- 
ered as an isotropic, homogeneous, compressible and viscous Newtonian fluid. The 
continuity equation is a manifestation of the fact that mass can neither be created nor 
destroved. The continuity equation states that the time variation of density within a 
control volume plus the mass entering and leaving the control volume is equal to zero. 
The differential form of the continuity equation for a compressible fluid and non-steady 


flow can be written as; 


—+V-(pl)=0. (1) 


For low speeds the density variation is small, therefore: 








and, 
Ve(pl’) = pVeV. 
Hence, for incompressible flow, the continuity equation can be written as; 


VeV=0. 


2. Derivation of the Navier-Stokes Equations 
For a compressible fluid, all primitive variables, density (p), pressure (p) and velocity 
(F) are functions of space and time. Newton’s Second Law states that the summation 
of all forces must be equal to the mass times the acceleration. 


- 


SF = Ma — (2) 
Considering an infinitesimally small fluid particle or control volume moving in a 
Cartesian Coordinate System, the right-hand side of equation (1) can be rewritten as: 


Ma =— (pV )dxdydz (3) 


D_ 
Di 


where 2 is the material derivative. 


Du Fic De De 


and J’ is the velocity vector, Which for a Cartesian Coordinate system is, 


Vu + yo twk (4) 


here u, v, ware the yelocity components along the coordinate axes. The external forces 
normally consist of the gravitational forces and the forces acting on the boundaries of 
the contro] volume, namely pressure and friction. All other body forces, such as 
electromagnetic forces will be ignored. For simplicity, the momentum equation only for 
the x-direction will be derived. while the derivation is analogous for the v and 2- 


directions. The x-component of equation (2) is: 








Max 2 (puldxdydz. (5) 


Figure (1) shows the normal and shear stresses on an infinitesimal contro] volume. 


Summation of the forces in the x-direction yields; 


a 


CT, Ct 
Surface Forces = o,dvdz + (ty, + ae dy )dxdz + (t,, + aor dz)dxdy 


—(o,+ = dx)dydz — t,,dxdy — t,,dxdz 


which reduces to, 


Co er . Ct. 
Surface Forces = (———* + ——— + — )dxdyaz. (6) 
éx cy éz 


Among the body forces only the gravity will be considered. Therefore, if f[ is the x- 


component of the gravity force then; 
Weight =f,(xy.z)p(xy,z)dxdvdz. (2) 


The sum of equations (6) and (7) are the external forces which are equal to the acceler- 
ation as stated by equation (5). After cancellation of the common term of volume 


(dx dy dz), the following force balance for the x-direction is obtained. 


COy + CT xy Cty 
Cx 





Se (PU) = Of + (- ) | (8) 


cy €2 


The next step is to express the stresses in terms of the primitive variables, i.e., velocities 


and pressure. First, the static pressure is defined as the mean of the normal stresses. 


pur (o,+o,+ o,) 


vf 


This equation can be algebraically rewritten as; 


o,=Pp +4 (20,- a, — 9;). 
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Figure 1. Normal and Shear Stress due to Friction 


In equation (9) the left-hand side is the normal stress at a point in the fluid. The first 
right-hand side term is the static pressure and the second right-hand side term is the 
deviation of the normal stress from the pressure due to viscous forces. Next a relation- 
ship between stress and rate of strain must be found. Isotropy implies that this relation 
between the components of stress and rate of strain is the same for every direction. The 
Newtonian fluid assumption means that this relationship is also linear. Referring to 
Figure (2) where o, and o, are resolved into diagonal components and equating the 


forces, the following force balance is obtained. 








ty (ay 2) - o,( = ) = a,( : ) =0 
ve \ 


; This equation can be rewritten as: 


, = y= (0x7 oy). (10) 












xy — plane 





Figure 2. Normal and Shear Stress 
A similar equation can be derived for the xz-plane. 
=+(6,-9; (11) 
Substitution of equations (10) and (11) into equation (9) results in; 


> 
On = PbS (t 2 — Tay): (12) 
The deformation of the initial shape of the fluid element (ABCD) to (A'BC’D) due to 
stresses in the y-direction is shown in Figure (3). The same figure also shows that the 
length of OA’ is as follows; 

(O.4') = Length = (u+ ay + H.O.T.)At. 


éx 
N 


The length change duc to stress is: 














Figure 3. Rate of Strain 


The strain is the change in length divided by the original length which is al. 2 . The 
strain rate is obtained by dividing this length by Ar. The same procedure can be re- 
peated for the y-direction to obtain the resultant rate of strain on a 45 degree plane due 
to o,. The shear stress on the 45° plane due to a, and a, is given by equation (13). A 
similar procedure for the zx-plane yields the analogous shear stress which is shown in 


equation (14). 











‘x =u SE ) (14) 


In the above equations, the proportionality constant yu, is defined as the coefficient of 
viscosity. Substituting equation (13) and (14) back into equation (12), and using Stokes’ 


Hypothesis; 
34 +2n=0 


yields equation (15) [Ref. 3: pp. 60-61]. 


Cu 2,éu , éy , Gv 
ee Wee ae tea) ie) 


In this equation the first term in parenthesis is the linear strain rate and the second term 
is the volumetric strain rate. To complete the derivation. the terms t,, and 1,, in 
equation (8) will be expressed in terms of the velocity components. From Figure (4) the 
rate of strain can be obtained [Ref. 19: p. 93]. 

The rate of strain on this element is _ Assuming that the variation of the 


rate of strain (;;) is small, the following expressions can be written. 


Cc) Cx 
.; Big ae eeeeoe 
Taking the limit of as Ar— 0, the rate of strain is given by: 
d, c Ay: 
a ee (16) 


Due to isotropy. t,, is equal to t,,.. Analogous procedures used to derive equation (16) 
can be repeated for the yz-plane and zx-plane, respectively, so that for a Newtonian 


Fluid equations (17) and (18) can be written. 

















Figure 4. Rate of Strain 


Cu Cy = 
ty =H +a) | (17) 
t= il > +—) (18) 


_ Finally, substituting equations (15), (17) and (18) back into equation (8), the momentum 


equation for the x-direction is obtained. 


c Cx Cx 
c Cu cy a cw Cu 
+— —=—+3——)|+— —+ > 19 
cy a cy Cx ] C2 [Hc Cx Cz )| At?) 





Similarly, the momentum equations for the y-direction and z-direction can be derived. 


These equations are; 


Dos) = pf, ee &| ye 29.7 
Dr (pv) = pf, ae alr Ee Ais 3 V ’)] 


cy ay oy 
+h [+H ]+2[u e+] (20) 
and, 
Plow = of, - 2+ [ue -25. r)] 


The unknowns in these last three equations are the primary variables; the density, the 
velocities and the pressure. (p,u.¥,4,p) . The momentum equations along with the 
continuity equation constitute the Navier-Stokes equations in the primitive variable 
formulation. The continuity equation for a Cartesian coordinate system is restated; 





é Cou py Cow 
o, fou, tor , ep 
Cr ex cy €2z 


=0. (22) 


Here the pressure is related to the density through the equation of state. 
p-pRT=0 (23) 


For an isothermal process and incompressible flow, equation (19) through equation (23) 
would be sufficient, but when temperature variations depend on density and pressure, 
the energy equation is also required. This is always the case for compressible flow where 
density depends on pressure and temperature. The energy equation expresses the bal- 
ance between heat and mechanical energy. The variation of viscosity due to temperature 
Variation may be obtained by an empirical viscosity law. The final result is a system of 
five partial differential equations with five unknowns; u, v, w, p, and p [Ref. 3, 19 } 
3. Derivation of the Energy Equation 

It is well known that energy can be neither created nor destroyed but it can only 

change in form. Therefore an energy balance exists for a fluid element in motion. This 
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energy balance is obtained through certain mechanisms which for a compressible fluid 
are determined by changes in heat content, total energy and mechanical work. Changes 
in heat content can be due to convection, conduction, friction and/or radiation. For the 
following derivation radiation is neglected because its effect is small at moderate tem- 
peratures. The energy balance for a contro! volume is expressed by the first law of 


thermodynamics. 
dQ _ dE, dW 
di at dt (24) 


First the variation of mechanical work or the contribution to work done by the external 
forces acting along the x-direction is derived. Again referring to Figure (1), the con- 
tribution to work done by each of the stress components is; 


eee ie cu £0. ee 
qi—6. = [ uo. + (u + “2 dx)(ox + ax) lavas 
which reduces to, 
ae. [ = (ua,) Jaxdvat. (25) 


Continuing with the same procedure for the other components of shear stresses; the y- 
direction and the z-direction, the total change in work due to normal and shear stresses 
can be written as shown in equation (26). 


a 


sae é ‘a 
-dl ae (uo, + VT. + wt,,) + ay (ut,, + voy + wry.) + ae (ut,, + vty, + vo.) |(26 


The total energy per unit mass within the control volume is the sum of the internal and 
kinetic energies, given by; 


72 


Total Energy =e + - ; (27) 


The variation in kinetic and internal energy for the control volume is shown in equations 
(28) and (29), respectively. 


dE internal = 4(pejaxdydz (28) 
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v2 
GE xinetic = dp a )dxdydz (29) 


Rewriting these two equations and summing them, the variation of total energy can be 


written as; 


dE 


v2 
dE = D- (pe + pA )dxdydz. (30) 


Changes in heat content due to conduction only are considered. According to Fourier’s 
Law the heat flux is proportional to the temperature gradient, so that the heat variation 
due to conduction can be written as; 


140 ~~ = aT 
Ad 4” er (31) 


Thus, by equating the amount of heat transferred into the volume with the amount of 


heat leaving the volume, the following relation is obtained; 


~k Lava: é (kD asyayds 


which gives the heat flux in the x-direction, 


= k at dxdydz. (32) 


Repeating similar procedures for the y-direction and the z-direction a final expression for 


the total heat variation is given by; 


dQ Oop OT seit é a 
7 -| = (k a) + eA ay \+ = (Se. dxdvdz. (33) 


Substituting equations (26), (30) and (33) back into equation (24), the general energy 
equation is obtained. [Ref. 3,19] 


D eg ly 7 sr LG Pi Dna Ne 
Fa a i aa a iD ar rr err 


A 


+— zx (uo, + Wty + WT,,) Boer ~ (uty xt to, + WTy») +>. E (ut, + vt, +wo,) (34) 
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B. CONSERVATION LAW FORMULATION 

The primitive variable formulation of the Navier-Stokes equations shown in the 
previous section can be put into conservation law form using vector identities. The 
conservation law form can be also derived by applying conservation principles on a 
control volume. Because physical insight is gained by this procedure the conservation 
law form derivation is outlined in the next section. 

This formulation stems from the fact that certain quantities (i.e. mass, momentum, 
energy) for a fluid in motion are conserved. Conservation implies that the flux of a 
quantity crossing a control surface and the net effect of internal sources results in a 
variation of the conserved quantity. These sources and fluxes depend on time and space 
as well as fluid motion. The fluxes are vectors for a scalar quantity and tensors for a 
vector quantity. Mass and energy are examples of a scalar quantity whereas momentum 
is a vector quanuty. Molecular motion and convective transport of a fluid contribute 
to flux. Molecular motion has the tendency to make the fluid homogeneous and has a 
diffusive effect. 

1. General Form of Conservation Law 

a. Scalar Conservation Law 
Considering a scalar quantity U within a control volume V’ , the time vari- 
auon of the quantity U is; 


“a 


—“/ UaV. 
Vv 


cl 


This should be equal to the incoming fluxes (F = ur) through a surface S (where n1 is 
the unit normal vector pointing outward), 


—| x +FaS=—J| F «ds 


plus any possible contribution from sources of U. 








Ht 


3h 


Figure 5. Flux Diagram 


The flux vector F has two components, a diffusive contribution and a convective con- 
tribution. The sources can be written as the addition of volume sources Q, and surface 


sources Os ; 
J Qual’ + J Os «ds 
so that the final conservation equation for the scalar quantity U is, 
é or Fe Aa osc 
= J eat =|OrdV +) OsedS—J F ods. (35) 
When using Gauss’s Theorem, equation (35) can be rewritten as; 


cU 


ve¢l 


dV4)VeFaV =| OV +) V6 OsaV 








For an arbitrary volume V the differential form of the conservation law is given by 
equation (36). 
cu V.(F-O 
a t ¥e(F -Qs)= Ov (36) 
b. Vector Conservation Law 

As stated earlier if the conserved quantity U is a vector, then the flux and 
surface source become tensors; F = Ue V , F, , and the volume source in turn becomes 
a vector, Q@,. An analogous derivation as in the conservation of a scalar quantity can 
be done for a vector quantity, whose integral and differential form are shown below. 


£ | ave[ FdS=| Ov +) O5-ds (37) 


él 


Cc _ => => _ 
“+ V6(F - Os) = Oy (38) 
In equation (37) the convective component of the flux tensor can be written in tensor 
form as: 

Fo, = yl} 
where r is the velocity vector. The diffusive component of the flux for a homogeneous 
system can be written as; 


Aa 


éu, 
Fp, =— pK on 

where « is the diffusivity constant. Equation (35) or (36) is the basic formulation of the 
conservation law for a general case. When continuity of flow properties is assumed (i.e. 
no shocks present), then equations (36) and (38) are valid. [Ref. 5: pp. 25-55] 
, 2. Equation of Mass Conservation 

In this particular instance the property U is mass and no diffusive flux is pres- 
ent, only convection. Therefore equation (35) can be directly written as; 


C —_ ~ 
<= | pa + | pV «dS =0 


or in differential form as in equation (39). [Ref. 5: p. 33] 
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(39) 


3. Equation of Momentum Conservation 
For this case the conserved quantity is momentum which is a vector. From 
Newton's Second Law it was mentioned that change of momentum is due to external 
volume forces and internal forces. Assuming a Newtonian fluid, the stresses can be 


written as; 
=> => 
o=—pl +t 


where / is the unit tensor, so that —p/ is the hydrodynamic pressure along the diagonal. 
The 7 term is the viscous shear stress tensor, equation (15), which is written as; 


ty = Hl(Gyy + 6x) — > (V+ F)6,). 


Referring to equation (37) and assuming that the external volume forces is zero the in- 
tegral form for the conservation of momentum is; 


> |> 


[pVaV + | pt «(F-ds)=[ 3 ods 
and applving Gauss’s Theorem, 
[prave) ¥-(r@navel sav 
Vet Vv V 


where © indicates the tensor product of two vectors. This can be written in differential 
form as shown below. [Ref. 5: pp. 40-50.] 
& (pi) +V«(oV QW’ + pl -7) =0 (40) 
4. Equation of Energy Conservation 
The quantity being conserved is energy, E and from the First Law of 
Thermodynamics the variation in energy must balance with the work of the forces acting 


on the system including any heat addition. The convective flux of energy can then be 


Written as; 


F-= pVE 
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where E is the sum of the internal energy plus kinetic energy. Definition of the diffusive 
flux term describes that diffusion of heat for a fluid at rest is due to molecular thermal 
conduction, and using Fourier’s Law of Heat Conduction, the diffusive flux term can 
be written as; 

Fp =kV+eVT 


where & is the thermal conductivity and 7 is the absolute temperature. Assuming no 
radiation, chemical reactions or work due to external forces, again the Q, volume source 
is zero. The net work done on the fluid by the internal shear stresses acting on the sur- 


face of the control volume is given by; 


Using equation (37) and substituting the quantities obtained above, the equation for 


energy conservation can be written as; 
2 pedv +) pEV doa [ATO dS 4] G+ ds 
= | ped + | pEV dS=| KTV aS +] (6 2V)0dS 
or in differential form as in equation (41). 


— (pE) +¥> »(V(Ep + p) — kVT—%-V) =0 (41) 


2|> 


Equation (41) can be rewritten as; 


De ene 2 CT é €T 
oe ee ee ae ee aoe + ££) 4 po 
PD, +P ax (Kay) t a = ay (KG tH 
where @® is called the dissipation function. Dissipation represents the heat equivalent 
of the rate at which the mechanical energy is lost during deformation of the medium due 
to viscosity. The dissjpation function is given by; 


o= 7-2 ty tse 4 Sey] (ee at 


Ox 


éu , cw . anny Cu , év , Cw % 
ey 363 Go GX See ey ee 





The system of partial differential equations given by equations (39), (40) and (41) can 


be written in compact vector notation as; 


é 


QI 


+V-0=0 (42) 


Y 


t 


where g is the vector of dependent conservative variables and QO is a vector composed 
of the nonlinear inviscid and viscous fluxes.[Ref. 5: pp. 45-50] 
5. Strong Conservation Form 
The strong conservation law form given by equations (39), (40) and (41) in 


vector notation can be written for a Cartesian coordinate system as; 


6g ay oo aioe aS a nae 
ai. 4 SE | GF, BG OR, eS. oT (43) 
cl Cx cy (oP Cx cy C2 


where q is the vector of conservative variables and E, F , and G are the flux vectors 


given by: 


p 
pu 
q=| pv (44) 
pw 
e 
pu pv pw 
pu’ + p pvu pwu 
E= pur F = pv’ +p Ge puy |. (45) 
puw pyw pw +p 
(p + e)u (pt+e)v (p + e)w 


The vectors of R, S , and 7 contain the viscous terms. When they are omitted, the Euler 


equations are recovered. 
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(T- V +4). 





The product term of 7 #’ is written in component form as below. 


(TV), = tut Tyt + TW 


(TeV), = tut tye + tw (47) 


(T+ V),=1,,u+ tay + 1,,W 


The heat flux vector g, is the heat transfer by conduction and can be written as; 


qe = ~kVT = —K(az, a, a3)" (48) 
where, 
ot I re 
aes) em Pr= 


In the above equations a denotes the speed of sound, Pr is the Prandtl number, c, is the 
specific heat at a constant pressure and e is the total energy per unit volume. [Ref. 20] 


Pressure and energy are related by the perfect gas law as follows. 
pHti- nfe -5W ++ »)| 


These equations can be transformed into different curvilinear coordinate systems in or- 
der to facilitate the numerical implementation. 

A coordinate mapping is introduced which allows the transformation of the 
equations of motion from a Cartesian coordinate, time varying. nonorthogonal coordi- 
nate system. The mapping is linked to the Cartesian coordinates as follows; 


& = Oxy 2.1) 
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n = n(x,y,2,1) 


f = C(x,y,z,!). 


The Cartesian coordinate system is the physical domain, and the transformed space is 
referred to as the computational domain. This computational domain is orthogonal with 
a uniform rectangular mesh so that unweighted differences can be taken to form the 
derivatives. 

The thin layer compressible Navier-Stokes equations are obtained from 
equation (43) by retaining the viscous terms only along the direction that is normal to 
the body. Also, the derivatives of the stress terms in the crossflow (i.e. y, z) directions 
are discarded. The thin laver formulation of the strong conservation law form of the 
governing equations for a curvilinear coordinate system (€,7,¢) along the axial, 


circumferential, and normal direction, respectively can be written as; 





Fy an - A N awe 
4, GAs (49) 
Cl cé cn CS Re ¢é% 
where g. F, G, Hand S are, 
‘ p pu 
pu pul + €p 
A 1 - 1 s 
q=7| pv Fawr) pyl + oyp 
pw pwl + €.p 
e (e+ p)U—Ep 
pV ply 
puV+nyp pull’ + Cp 
r ] ; an 1 . 
G TT pxV+nyp H= 7 pri’ + Cyp 
pwV +p pwlhl’ + Cp 
(e+ p)V—np (e+p)— op 
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0 
. pmyu; + (u/3)m0 
S= > pm Vv, + (u/3)m 6, 
pm, w, + (u/3)m2C, 
umym, + (u/3)m + 2(0,u + Cv + 0,w) 


Furthermore, it is defined that; 
2. yd 62 
m= (4a + by + - 
my = Cu + Cyr + Uwe 


Pe | 
inf Pe se K , Ga 
m,=(u'+v + w)/2 + >> ("ar ) 


and U, V, and WV’ are the contravafiant velocity components given by, 
Usut,+ve,t we, te, 
Veun, ty, + wy, t ny 
Weul,+wytwt,t,. 


Again analogous to the previous derivations the pressure is related to density and total 


energy through the equation of state for an ideal gas. 


C. NUMERICAL IMPLEMENTATION 
1. The Numerical Algorithm 

The solutions over a strake-delta wing configuration resembling a modern 
fighter aircraft planform will be presented in the last part of this thesis. Even though the 
main effort of this work was not the numerical solution of the governing equations (i.e. 
the compressible Navier-Stokes equations), the technique used for the numerical imple- 
mentation is briefly described in the following paragraphs. 

The numerical scheme used for the solution of the governing equations is based 
on 4a finite difference discretization of the thin layer Navier-Stokes equation [Ref. 6]. 
The numerical integration was performed using a partially flux-split numerical scheme. 
Upwinding was performed in the main flow direction using flux vector splitting while 
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central differencing was used in the other two directions. The factored form of the re- 
sulting algorithm is as follows; 


[1+ hd3(A*)" + hd-C” — hRe~"5, I'M" — Dy | 


x [1+ hdf(A7)’ + A5,,B" — Dy, Aq” = : 
— A(32[ (F*Y" — FL] + LF)" — FL] + 5,(G" — Gq.) + ¢(Hy — Hon) + > 8(S” — 5.0) 


= Aq" _ Yoo) 


The explicit dissipation D, was used along the directions where central differencing was 
applied. The implicit dissipation term D, was added for numerical stability. Steady 
aerodynamic flows at subsonic flows (M = 0.2) do not contain shock waves and can 





be quite well predicted by a central difference scheme that is augmented by these dissi- 
pation terms.[Ref. 20 } 
2. Turbulence Model 

Simulation of high Reynolds number flows is obtained by the solution of the 
Reynolds averaged Navier-Stokes equations. These equations have extra unknowns and 
are commonly called the Reynolds stresses [Ref. 3}. The relations between the Reynolds 
stresses and the mean flow quantities is the well known closure problem. In practice 
some turbulence model is used which relates the Reynolds stresses with the mean flow 
quantities. The turbulence model selected for this research was an algebraic eddy 
viscosity. This model is the Baldwin-Lomax model as modified by Degani and Schiff to 
treat three-dimensional separated flows [Ref. 21,22]. 

The turbulence is simulated in terms of an eddy viscosity coefficient u,. The 


coefficient of viscosity and the heat flux term in the Navier-Stokes equations are re- 
H 

R Pr , 
developed by Cebeci with modifications that allow for the locating of the boundary layer 


placed with » + yu, and > + respectively. The turbulence model is similar to one 
{Ref 23}. A two layer algebraic eddy viscosity model is used where the Prandtl-Van 
Driest formulation is used in the inner region and the Clauser formulation is used in the 
outer region [Ref. 21]. The inner region is any normal distance from the wall, y , that 
is less than or equal to J, + If this is the case then y, is defined by the following ex- 


pression; 


(H,)inner = pl’ | @ | \ 
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where, 





—)* 
l= wf — exp( a } 


and, 


ev vt cw 2 cw Gu 
Se ea de ae eee We ee mer 
Vv éy  éx ren ez cy y+ 





If y is greater than },,,.,.-, then u, is defined by; 


(Udouter = KCopF wareF ries?) 
where A is the Clauser constant, C,, is an additional constant, and for boundary iayers, 


Fyake = Ymaxl 


max 


or for wakes and separated boundary layers, 


Uni 
Di 
Fare = Cwr)'max Fouax 


In the above equation U’,,, is the difference between the maximum velocity at Jn. and 


the minimum velocity in the profile. The quantities of y,,,, and F,,,, are calculated using 


. =" 
Fo) = stot] 1 esp 7 )} 





The function Fy,.() is the Klebanoff intermittency factor and is defined as; 


Crten¥ z= 
Frias) = [ + 5.5( eee ‘| : 


J max 


All other parameters are constants determined empirically and given in [Ref 21]. The 


use of this model eliminates the need for finding the edge of the boundary layer and re- 
duces one of the sources of error in the Navier-Stokes solutions. 
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III. SURFACE GRID GENERATION 


The focal point of this work was the generation of the computational mesh over a 
three-dimensional strake-delta wing configuration that models a modern fighter aircraft 
planform. Therefore, in the next two chapters the surface and field grid generation 
procedures are described. The considerations which must be taken into account in order 
to construct a surface grid suitable for the subsequent generation of the field are dis- 
cussed. Finally, the various field grid generation methods are discussed and two different 
approaches which were used to generate the field grid over the strake-delta wing con- 
figuration are described. 

The first step to be taken in establishing a finite difference or finite element scheme 
for solving a system of partial differential equations is to replace the continuous domain 
by a finite mesh, commonly known as a grid. Grid generation is one of the central 
problems in the procedure to obtain a numerical solution. A well constructed grid 
greatly facilitates the numerical solution of a system of P.D.E.s. On the other hand, an 
improper grid choice may lead to instabilities, inaccuracies and’or lack of convergence. 
Numerical grid generation is a procedure for the orderly distribution of observers over 
the physical field domain in such a way that efficient communication among the ob- 
servers 1s possible. Also, it assures that all physical phenomena of interest in the entire 
ficld may be represented with sufficient accuracy by this finite collection of observers. 

Grid generation for two-dimensional domains is relatively simple and may be 
achieved with purely algebraic techniques, even for relatively complex domains [Ref. 
24). In addition, for suitable geometries of the boundaries conformal mapping techniques 
may be used. Conformal! mapping techniques have the advantage that they are relatively 
simple and inexpensive [Ref. 25 : pp. 488-490][Ref. 4: pp. 7-56]. They also preserve grid 
orthogonality, but their use is limited to domains with simple boundaries where a con- 
formal transformation between the physical domain and a simpler transformed domain 
may be readilv defined. 

The generation of a computational grid for a three-dimensional domain, however, 
presents greater difficulties. For a limited class of external and internal domains it is 
sometimes possible to fill the entire three dimensional domain with a sequence of two- 
dimensional plane grids that will constitute the entire three-dimensional field grid. An 
application of this idea is shown later for the construction of the field grid over the 
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double-delta wing. In many instances however, it is either difficult to decompose the 
three-dimensional domain into a sequence of two-dimensional domains, or it is prefera- 
ble to construct a purely three-dimensional mesh. Of course, both the complexity and 
computing time of a three-dimensional grid generation method will be higher. In any 
case, the definition of the surface boundaries must be done precisely and accurately. 

Before any work can be started on a three-dimensional field grid over a body one 
must first define the surface geometry of the body. The definition of the body’s surface 
and its quality is imperative to the success of the field grid. There exist many avenues 
to create the surface grid, some include algebraic techniques, cubic Hermite functions, 
Bezier curves or Non-Uniform Rational B-splines (NURBS) [Ref. 4: pp. 237-249] [Ref. 
25: pp. 497-503]. Each of these techniques has its advantages and the exact method that 
will best fit a particular surface will vary. The availability of accurate data for de- 
scription of the surface geometry will also play a major role in the generation of the 
surface grid. If all surfaces can be defined in terms of equations, then an algebraic 
technique might prove to be the most efficient. Whereas, if the surface is very complex, 
as is the case for actual aircraft surfaces, all that is available are two-dimensional cross- 
sections and a curve fitting technique will have to be used. Whatever the method, it is 
of utmost importance that the surface grid generation program be written in such a way 
that it will enable maximum flexibility in the number of grid points and their distrib- 
ution. This early concern and respect for versatility will pay large dividends upon sub- 
sequent generation of the field grid. For even after the surface grid has been completed, 
an interactive trial and error process of changing the surface grid will be required to 
achieve an effective field grid. 


A. DOUBLE-DELTA WING SURFACE GRID 

The dimensions of the double-delta wing model are shown in Figure (7). From this 
figure it can be easily seen that most of the surfaces can be defined by linear relation- 
ships with the only exception being the NACA 64-005 cross-section. For this reason 
algebraic grid generation on the surface was chosen. For a linear relationship and ana- 
Ivtically defined points no advantage is gained by using a curve fitting method. The part 
of the wing that contains the NACA 64-005 airfoil cross-section required special treat- 
ment. Generation of the surface grid over the airfoil cross-section would require a curve 
fitting technique. Much work has been done in this area and NURBS can produce ex- 
cellent results in approximating airfoils. The main advantage of this technique is that 
it provides the flexibility of modifving the cross-section or shape of the surface by simple 
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changes of user specified parameters. The disadvantage is that the complexity is higher 
and the redistribution of the points that represent the airfoil is also more difficult. Be- 
cause the cross-section is a NACA airfoil whose contour shape can be well approxi- 
mated by straight lines, the use of a purely algebraic technique for the entire 
double-delta wing surface grid was utilized. 

Once this decision was made, the areas containing singularities had to be identified 
and the program for the algebraic grid generation had to be written. Because the surface 
grid is used as initial or boundary conditions for the generation of the field grid, much 
care had to be taken to avoid any singularities that would propagate into the field grid. 
In addition, special care must be taken at the regions of sharp corners such as the lead- 
ing edge. In these areas it is not possible to maintain field grid orthogonality and the 
location of these acute angles can be seen in Figure (7). In general, these areas are found 
on the entire leading edge, at the apex and at the rectangular edge near the wingtip. 
These corners had to be approximated by “rounding” off these areas with a radius that 
was very small. The radius used was 0.001% of the chord length, so to the naked eve 
the surface grid appears to be a sharp corner. This rounding of acute angles allows the 
field grid to maintain orthogonality which is a desirable feature for subsequent numerical 
implementation; see Figures (11), (14), (17), (20) and (22). Also, high grid resolution is 
provided at the same time in these areas where the change of the flow field variables is 
expected to be rapid. The methodology for generating the source code that would 
compute the grid points was to progress from the nose in an axial direction through to 
the wake. The grid points were essentially generated for a two-dimensional cross-section 
in the yz-plane, then an incremental step in the x-direction was made and again the grid 
points for a new yz-plane were computed. The source code was written in five logical 
sections that defined regions of the wing with similar cross-sections. These five sections 
were the apex, the strake, the wing, the trailing edge rectangular section and the wake. 

1. The Apex 

Special care had to be taken in the modeling of the nose region. The apex of 
the wing is a single point that transitions to a diamond shape cross-section. Taking into 
account that smoothness has to be maintained, a hemisphere mav be used to provide 
smooth transition between the singular point of the apex and the diamond cross-sections 
at the nose, see Figure (11). The radius of this hemisphere is 0.001% of the chord length 
and allows for a smooth transition. The radius of the sphere, the number of grid points 
for the axial (x-direction) and the circumferential (v-direction) were inputted by the user. 
An incremental angle was determined for both the xz-planes and yz-planes which then 
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enabled the computation of the grid points. The yz-plane cross-sectional grid points 
were then calculated using incremental yz-plane angles as the x-location progressed 
downstream using the incremental xz-plane angles. 

2. The Strake 

The main concern in defining this section was that the leading edge has a corner 
that is relatively sharp and would preclude a field grid that is orthogonal. Therefore, the 
sharp leading edge was approximately rounded as shown in Figure (12) through Figure 
(17). The computations here involved a user specified radius that was the same as the 
one in the approximation of the apex. This radius was maintained to allow a smooth 
transition from the sphere of the apex to the diamond cross-section of the strake. The 
surface grid generator code provides the versatility to change the number of grid lines in 
this radius which is kept constant for every cross-section. This issue becomes important 
as the ratio of radius to wing span drastically changes between the apex and junction 
of the strake and wing. Again similar logic to the one used to define the apex was used 
here. The difference being that the increment of the x-coordinate was computed de- 
pending on the number of grid lines along the x-direction used to define the strake part 
of the body. Simple relations from analytic geometry are used to determine the y and 
z-coordinate as a function of the x-location. 

The distribution or clustering of the grid lines will be discussed in more detail 
later in this section. It is important to mention that the distances between successive 
grid points, along the x and y-directicns was determined by calling a subroutine. This 
allows to experiment with many different distributions with a simple change of input 
parameters. 

3. The Wing : 

For the wing section three areas required special attention. The first was like 
the strake, in that the leading edge forms a corner which unlike the strake was not as 
sharp. This was because of the NACA 64-005 cross-section has a finite curvature at the 
leading edge. Part of the leading edge did not require rounding and the number of grid 
lines approximating the edge was reduced, see Figure (15) through Figure (17). The 
NACA 64-005 cross-section is generated by a subroutine which requires as input only 
the normalized root chord length of the airfoil, x, = x,(). The area of the wing spanning 
between the wing centerline and the part of the wing having a NACA 64-005 cross- 
section was defined by linear interpolation. Some sort of curve fitting method could 
have been used, but since the wing is.thin and the distance is short, a linear approxi- 
mation was assumed to be sufficiently accurate. 
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4. The Trailing Edge Rectangular Section 

This particular section required the most challenging surface grid definition in 
the early stages. This was primarily because the wingtip has a variable finite thickness 
that depends on the x-location, see Figure (18) through Figure (20). The actual calcu- 
lation of the grid point locations was simple, but the number of grid lines at the wingtip 
region had to change due to this variable thickness. This required to take grid lines out 
of the edge and redistribute them back onto the upper and lower surfaces. This is the 
reason why the grid lines in this section, when viewed in the xy-plane appear staggered, 
see Figure (9). The trailing edge has a finite thickness (0.002% of chord length) and thus 
avoids any singularities that unnecessarily complicate the subsequent numerical imple- 
mentation and the generation of the field grid. 

5. The Wake 

For the numerical implementation, an extension of the far end of the computa- 
tional domain of 2.0 - 3.0 root chords beyond the body is required. The wake was rela- 
tively simple to generate because all that varied was the x-coordinate. All the yz-planes 
remain constant from the trailing edge to the end of the grid. Examples of this cross- 
section can be seen in Figure (21) and Figure (22). The wake extends for 2.0 root chord 
lengths bevond the wing trailing edge. This length was selected because it allowed for 
a smooth transition from the wing to the wake for a given number of x-direction grid 
lines. 


B. DISTRIBUTION PARAMETERS 

Flexibility in the distribution of the surface grid points in both the x-direction and 
the y-direction, which are shown in Figure (10), is important for the surface grid. The 
- foremost problem is to ensure that the distance between successive grid points makes a 
smooth transition. Of course, the spacing between the surface grid points could have 
been made the same, but this fs impractical because the total number of grid lines would 
be excessive due to the small radius that was used to approximate the corners. Therefore 
a distribution of grid points must be developed that is very dense at the corners and 
sparser in the other regions. High grid clustering is also required in areas where steep 
gradients in the flow-field are expected, such as the leading edge where the leading edge 
vortices appear. The distribution in the y-direction can be seen in Figure (9) and Figure 
(10) for various cross-sections of the wing and the distribution for the x-direction can 
be seen in Figure (10). The y-direction has a high grid clustering around the leading edge 


which becomes sparser near the centerline of the wing. This was the general procedure 
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followed for the y-direction distributions in all cross-sections. The distribution along the 
x-direction required a higher density in the nose region and sparser distribution at the 
area near the end of the wake. A subtle change to a higher density occurs where the 
wing has a geometry change as can be seen in Figure (9). 

Stretching of the grid points along a coordinate direction can be obtained by using 
simple algebraic functions such as linear, exponential mapping or trigonometric func- 
tions. Use of these functions allows for a smooth transition from sparse grid densities 
to high grid densities. A quadratic function was first attempted but this resulted in a 
distribution that became less dense too fast and would spread the grid points out to an 
excessive amount near the centerline. Next a linear stretching function was used and this 
gave much better results but did not allow a “smooth” transition from the high grid 
density region to the region with sparser grids. This effect was more pronounced along 
the y-direction. The linear function allows a more constant distribution over the whole 
Wing in the x-direction. The linear equation shown below was used for the linear 
stretching: 


Xj = CX; 


where the user specifies the parameter c depending on the desired degree of stretching. 
A value ofc = 1.0 would result in an equidistance spacing, whereas a c = 2.0 would result 
in a high clustering of grid points near one or both ends. Difficulties were not en- 
countered in the x-direction because the transition from the low to high density distrib- 
utions Were not as extreme as in the y-direction. By changing the linear stretching 
parameter, a smooth transition from more to less dense areas was achieved. The grid 
stretching in the axial direction required different values of c depending on the wing 
section; for example, c = 1.005 was used prior to the trailing edge and c = 1.205 after the 
trailing edge of the wing. The effect of these constants can be seen in Figure (9). Each 
representative cross-section had to be investigated and a constant assigned that allowed 
a smooth transition from one section to another. These constants were determined by 
a tria] and error procedure, but experience gained by many iterations expedited the 
process. From the many iterations for the surface grid alone, an appreciation for the 
flexibility of the source code was gained. 

To resolve the y-direction distribution, the stretching function first attempted was 
exponential which proved to be inadequate. The exponential stretching is obtained by 


using the following expression; 








Set = Vi 


where c and s are parameters chosen by the user to produce a desired distribution. The 
dimensions of the radius used to approximate corners are so small compared to the 
characteristic dimension (i.e. the chord length) that for the desired number of grid points 
in the y-direction the transition did not occur smoothly. Finally, a sinusoidal distrib- 
ution produced better results. The equation below was used to obtain this distribution; 


Ye=csinb 


here c is a user specified parameter and b is allowed to incrementally change over a 
specified range of angles in order to obtain the desired section of the sine curve. The 
distribution produced the best results for a sine curve segment from 0 to 45 degrees. 
Freedom was written into the source code to use constants to finely adjust the distrib- 
ution but were not required because of sufficient results without them. 


C. PROGRAM FEATURES FOR THE SURFACE GRID 
The source code for the surface grid can be seen in Appendix E. The main concern 

in the construction of the source code for this problem was to give the author maximum 
flexibility in the generation of this surface grid. Listed below are some of the features 
that can be easily changed via an input file. 

¢ The number of grid points in the x-direction at five different sections. 

e The number of total grid points in the y-direction. 

e The radius used in the approximations of sharp edges. 

e The number of grid points used in the radius for the corner approximations. 

e The sweep angles of the strake and the wing. 

e Maximum widths of the strake and wing. 

e Lengths of the strake, the wing, the rectangular section and the wake. 

¢ Distribution constants at five locations in the x-direction and at six locations in the 

y-direction. 

This flexibility in the surface grid paid a major dividend in the subsequent generation of 
the field grid. This is because the surface and field grid generation is an interactive 
process that usually requires changes in the surface grid. Another feature is that the 
distribution functions are written as subroutines. This allowed the author the flexibility 
of trying different functions based on the geometry and desired gradients of distribution. 
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Another subroutine is the calculation of the grid points that are part of the NACA 
64-005 airfoil cross-section. This enabled the changing of this cross-section by merely 
changing two lines of the data statement. The subroutines that are included in Appendix 
E are linear, quadratic, ex.onential and sinusoidal. Again it is emphasized that the 
source code for the surface grid generation must be as versatile as possible. 

This program can be used to generate a surface grid for a wing with dimensions that 
are different from the one chosen by the author. But because great care needs to be 
taken in grid point distribution a change in the dimensions would most definitely require 
an adjustment of radius, the number of points, distribution constants and even distrib- 
ution functions. A close examination of every grid constructed, either visually or com- 
putationally, must be completed to ensure that no intersection of the grid lines occurs. 
The program that was originally written was continually revised to permit an adequate 


construction of the field grids. Included in Appendix F is the final program that was 


utilized for the generation of the surface grid for the final spherical field grid topology. 








IV. FIELD GRID GENERATION 


In this chapter different numerical techniques for the field grid generation are pre- 
sented. Their advantages and disadvantages as far as grid generation and numerical 
implementation are explained. Most available field generation techniques require as in- 
put a surface grid which must be constructed before the field grid generation. The gen- 
eration of the field grid is directly dependent on the distribution of the grid points on the 
body surface. There are several ways to generate field grids, a few of which include 
methods involving solutions of Elliptic Partial Differential Equations (P.D.E.), Parabolic 
P.D.E.s, or Hyperbolic P.D.E.s. For a limited class of problems algebraic methods can 
also be used. In the following paragraphs the various grid generation methods based 
on the solutions of P.D.E.s will be described. The hyperbolic method, which was used 
to generate the field grid over the double-delta wing will be explained in detail, whereas 
only a brief description of elliptic and parabolic techniques will be given. 

The classification of P.D.E.s into hyperbolic, parabolic or elliptic type is obtained 
fro1a the general form of the quasi-linear second order P.D.E., given by: 


Quy, + bu,, + Cuyy + du, + eu, + fu=g (1) 


here u = u(xy') is the dependent variable and the coefficients a,b,c,d,e,f and g are func- 
tions of x andy. The type of equation (1) is determined from the sign of the quantity 
b? — dac as follows. 


b?~4ac<0 (Hyperbolic) | (2) 
b? —4ac=0 (Parabolic) (3) 
b?—4ac>0 (Elliptic) (4) 


Each type of equation, hyperbolic, parabolic or elliptic has certain characteristic prop- 
erties which can be successfully utilized for the grid generation in two and three- 
dimensional domains [Ref. 4: pp. 188-277]. For example, the solution of an elliptic 
equation in the interior of a domain depends on the specification of data over the entire 
boundary. Therefore, when a grid is generated by the solution of an elliptic P.D.E. all 


the boundary data must be specified. 

















The main feature of a hyperbolic type of problem is that the solution starts from an 
initial condition and propagates in time along certain directions known as the charac- 
teristic directions. Utilization of this property allows the construction of grids over 
contoured lines or surfaces by propagating in space the initial information provided by 
these lines or surfaces. The grid generation method must be carefully chosen to facilitate 
the type of grid desired. Each type of grid generation method will require some addi- 
tional data to allow for a suitable solution. To some extent this additional data may 
determine which type of grid method should be utilized. If x and y are spatial coordi- 
nates (which is true for the present case) then the additional data will be the boundary 
conditions. If x and y represent time then the additional data will represent initial con- 
ditions . 


A. ELLIPTIC GRID GENERATION 

Field grid generation obtained by the solution of an elliptic set of equations requires 
specification of each and every boundary point of the closed domain where the grid will 
be generated. The inner boundary is simply the body surface grid and the outer 
boundary 1s a user specified shape. The body surface must be specified exactly. How- 
ever, there is some flexibility in choosing the shape of the outer boundary. Another re- 
quirement of this method is that the curvilinear coordinates must be constant or 


monotonic on the boundaries. If any extrema of the curvilinear coordinates exist in the ~ 


interior of the physical region then overlapping of the grid lines will occur. When using 
an elliptic method, initial boundary slope discontinuities are not propagated into the 
field. This feature of elliptic grid generators tends to make the grid very smooth. The 
large computational time requirement for the solution of the elliptic system of P.D.E.s 
can be a disadvantage. The simplest form of an elliptic P.D.E. is the Laplace equation; 


Vei=0 (5) 


where ;= 1,2 for two-dimensional grid generation and i= 1,2,3 for three-dimensional 
‘grid generation. The effect of the Laplace operator is that a very smooth grid is 
produced which becomes equally spaced away from the boundary. The Laplacian also 
guarantees one to one mapping of the coordinate system. This method will have the 
effect of making the grid lines more closely spaced over concave boundaries and sparse 
over convex boundaries.[Ref. 4: pp. 188-228}[Ref. 25: pp. 503-510] 

Another approach to generate the field grid is to solve a Poisson system of 
equations. This system has the following general form; 
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Ve! = pl. (6) 


The forcing term P' can be used to control the spacing and orientation of the grid lines. 
This control can be extended to move the intersecticn slope of the grid line with the 
boundary. When P'- 0 the grid lines tend to become equally spaced, i.e. to approach 
a grid obtained from the solution of Laplace’s equation. The forcing term P’ can also 
be used to enhance grid orthogonality. Orthogonality of the grid lines close to the body 
Surface grid does not occur normally in elliptical solutions. The main advantage of a 
Poisson type grid generator is that orthogonality control of the grid lines can be main- 
tained at the expense of complex, lengthy and expensive calculations. There are several 
elliptic grid generators in use today that maintain grid orthogonality.[Ref. 4: pp. 
193-236.] [Ref. 26] 


B. HYPERBOLIC GRID 

The hyperbolic method involves marching in space in a time-like fashion of the 
boundary information, i.e., a surface grid. This method is suitable for external flow 
problems where the exact location and shape of the outer boundary is of no vital im- 
portance. One major advantage is that computationally this method is efficient; in ad- 
dition, orthogonality of the field grid is preserved. Hyperbolic methods are usually one 
to two orders of magnitude faster than the elliptic methods because of their noniterative 
time-like marching nature. Control of the grid lines is somewhat restrictive but specifi- 
cation of the cell volume can result in the avoidance of overlapping grid lines, especially 
in concave areas. Overlapping of the grid lines is not allowed because singularities are 
propagated into the field, so great care must be taken to avoid these when constructing 
the surface grid. Because the characteristics of the hyperbolic method include 
orthogonality preservation and computational efficiency, a hyperbolic grid generation 
method was selected as opposed to an elliptic grid generation method. Great care was 
taken in the generation of the surface grid to remove any singularities such as sharp 
corners, because rapid transitions on the surface geometry usually produce intersecting 
grid lines of the field grid. The wing configuration did not contain any severe concave 
surfaces which would cause intersecting of the grid lines. Two different grid topologies 
were examined for the grid generation of the field grid over a double-delta wing. i.e., 


evlindrical and spherical. First the cylindrical grid generation procedure is presented. 
[Ref. 4: pp. 272-276}[Ref. 25: p. 503] 














1. Cylindrical Grid Generation 

The cylindrical type grid produces an H-O configuration shown in Figure (23). 
The grid in this figure is called H-O because when the wing is viewed from above or from 
the side the grid appears to have a “H” shape. When the grid is viewed from a nose-on 
viewpoint (i.e. the vz-plane), then the grid has an “O” shape, see Figure (24). It is usu- 
ally easier to generate a cylindrical grid than generating a spherical grid because the 
three-dimensional cylindrical grid is an assembly of planar two-dimensional sections. 
This combination of two-dimensional planes starts at the nose and progresses aft to- 
wards the trailing edge, see Figure (25) and Figure (26). Various methods may be used 
to generate these two-dimensional grids on the planar cross-sections. Here, a two- 
dimensional hyperbolic grid generator was used to generate the plane O-type grids at 
various locations along the axial direction. The cylindrical grid has a singular point at 
the apex, Figure (25), where special care must be taken during numerical solution for the 
computation of the transformation metrics and the application of boundary conditions. 
Along the singular line starting from the apex and extending upstream to the beginning 
of the domain all the grid points collapse on to a single point. 

Before the surface grid data points generated for a generic surface can be used 
for a cylindrical grid, a modification was required at the nose of the grid. As can be seen 
in Figure (7) a singularity is present at the very first point at the apex. To avoid the 
problems that this point can cause, the singular point at the apex was omitted. The field 
grid was then generated to have an annular type appearance in the yz-plane cross- 
section, see Figure (29). Because the radius is so small, the effects due to the inaccuracy 
of its surface definition are negligible. 

The process of improving the quality of the field grid was completed using the 
computer graphics program PLOT3D [Ref. 27]. This graphics package was designed to 
facilitate visualization of the field and surface grids and flow fields of computed and ex- 
perimental results. This same program was utilized to correct the surface grid during its 
development process. Because of the early concern regarding the surface grid and 
knowledge of hyperbolic grid shortcomings, only minor adjustments were required in the 
surface grid. These adjustments required a redistribution of surface grid points in the 
spanwise direction. 

The initial surface grid had a linear distribution that was deemed inadequate just 
by visual inspection and from physical considerations of the flow field character at the 
leading edge region. A quadratic distribution was subsequently used which appeared to 
yield a better distribution of the grid points. Not until the field grid was generated was 
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it clear that this distribution was still insufficient. The necessary changes in the field grid 
were to smooth out the distribution near the leading edges and to increase the grid 
density where vortices were most likely to occur. Finally, a sinusoidal distribution along 
the spanwise direction gave the best results. Examples of various cross-sections of the 
surface grid distribution for the wing can be seen in Figures (29), (31), (34) and (37). 
At the trailing edge the same reasoning as for the nose was carried out and the grid was 
not collapsed onto a single line but instead retained a finite thickness. Asa result a very 
small but finite thickness wake was generated. 

The field grid was completed by extending the grid from the wing apex to a lo- 
cation 2.0 - 3.0 chord lengths upstream, where the freestream conditions can be applied. 
The reason for this ‘s that the flow field is affected upstream by the presence of the wing. 
This addition to the grid was completed by repeating the very first annular shaped grid 
and by changing only the x-location. The yz-plane remained the same and it was re- 
peated as the x-distribution gradually increased its Ax spacing until it reached a user 
specified x-location upstream. The short program for this additional grid can be seen in 
Appendix F and examples of the entire field grid (cylindrical type) can be seen in Figure 
(24) through Figure (41), Appendix B. The completed cylindrical field grid dimensions 
are 110x240x68. 

The purpose of the grid generation is to facilitate a numerical solution to a 
system of P.D.E.s, and for accurate solutions to occur certain areas require special 
treatment. A problem that occurs in the computation of derivatives in the apex region 
is that differences are taken between points that may have different flow characteristics; 
1.€. One point may be in the boundary laver and the next point may be outside of this 
regime. This is the reason why clustering of grid points near the apex is required to 
avoid as much as possible these inaccuracies. Clustering of the grid points normal to the 
surface of the wing is naturally required to resolve the velocity gradients in the viscous 
boundary layer. For the most part though, the grid is aligned with the main flow direc- 
tion and a numerical scheme that uses upwinding can produce accurate results.[Ref. 20] 

2. Spherical Grid 

This particular method of grid generation produces a C-O type grid as can be 
seen in Figure (24). This method is more complex than the cylindrical one because the 
entire three-dimensional grid must be generated simultaneously. The spherical grid 
topology has one singularity that is located at the apex and propagates upstream, see 
Figure (42) and Figure (43). The spherical grid is also aligned with the main flow; 
therefore, an upwinding scheme can be used for flow field solutions. One disadvantage 
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of the spherical grid topology is that the visualization of the flow is a little more difficult 
than for the cylindrical grid because the € = constant grid surfaces are not planes but 
three-dimensional surfaces. The flow field solution tends to converge faster on a spher- 
ical grid than in the cylindrical grid. 

The hyperbolic grid generator for this work utilized a cell volume hyperbolic grid 
generation scheme. A coordinate transformation to the computational domain (€, ”, f) 
was performed where the body surface was the boundary condition ¢(x,y,z)=0. The 
field grid is obtained by the solution of the following nonlinear system of P.D.E.s; 


XX +I + 2y2 = 0 
XV yZe Tt AXV 32 XV Zs — XV eZy — Xyy2e — XVyZs = AV. 


where the initial condition at { =0 are the x,y , and z coordinates of the body surface 
{Ref. 20]. The first two equations are the relations that preserve orthogonality with re- 
spect to an outward normal vector ¢. The third equation is a user specified volume 
parameter that controls the cell size and normal spacing of the grid points. The grid is 
generated by “marching” in the { direction and the system of P.D.E.s is solved by an 
approximate-factorization, noniterative, implicit finite difference scheme. Even though 
grid orthogonality and smoothness are maintained the quality of the field grid is quite 
sensitive to the quality of the surface grid. Control of grid clustering along the normal 
to the surface direction is provided, but there is no accurate control in the location of 
the outer boundary due to the marching type solution. The outer boundary for the 
purposes of the present work is not crucial as long as it extends 2.0 - 2.5 root chord 
lengths away from the body. 

The three-dimensional hyperbolic grid generation is very sensitive to the surface 
grid definition, since the surface grid distribution is propagated in space to generate the 
three-dimensional mesh. This sensitivity to the initial conditions is the reason why the 
spherical grid generation presented more difficulties than the cylindrical grid generation. 
The apex singularity in conjunction with the sharp angles created most of the problems 
in the grid generation as far as preservation of orthogonality is concerned. Because of 
this point singularity, a blunting of the nose was first attempted. This nose region also 
resulted in a reduction of grid lines in the x-direction to 110. For the C-O configuration 

the first grid line extends upstream and therefore it is not necessary to add axial grid lines 
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as in the cylindrical grid case. For effective field grid generation, distributions in all three 
directions had to be adjusted near the nose. After many iterations a solution was gen- 
erated that required only minor adjustments. These adjustments were made by writing 
a small program that would algebraically adjust the x-distribution at the nose singularity. 
Algebraically fixed were the second and third grid planes in the x-direction. This tech- 
nique is sometimes the only option available when the grid requires minor adjustments 
and the changing of usual parameters produces negative results. Although a suitable 
grid was constructed, an alternative method of allowing the apex to collapse to a sharp 
point was attempted in hopes of an even more suitable field grid. Also at this time sol- 
utions were being generated for the cylindrical grid and it was observed that the field grid 
required a higher grid density in the wing area to facilitate better definition of the vortex 
that develops over it. For this reason the number of grid lines in the axial direction was 
increased to 160. By allowing the apex to converge linearly to a point resulted in a 
three-dimensional mesh that had been unsurpassed up to this point. The disadvantage 
of allowing the strake to linearly collapse to a point was the original y-direction 
sinusoidal distribution (0° to 45°) of the surface grid was now insufficient. A solution 
to this problem was to allow the y-direction distribution near the apex to have a 45° to 
60° distribution. This distribution results in a spacing of the grid lines that is early 
linear. Then the distribution was allowed to change incrementally as the x-location 
changed to achieve a 45° to 90° sinusoidal distribution at the strake and wing junction. 
This change was only required in the y-direction. The conscious decision to allow a 
singularity at the apex did not affect the quality of the field grid because the topology 
of the spherical C-O type grid requires that this area collapse to a singular line.. 

Other areas of the wing did not require adjustments because the surface grid did 
not propagate any problems into the field grid. Cross-sections of the grid can be seen 
in Figure (45) through Figure (56). The only other adjustment made was to stretch the 
entire grid to a suitable distance from the wing. This also was done with a small pro- 
gram that operated on the data file that was generated from the hyperbolic grid genera- 
tor. While the final grid appears more uniform throughout the entire field, slight 
deviations in the orthogonality to the surface did occur. The final grid dimensions for the 
spherical field grid were 160x240x68 and can be seen in Figure (42) through Figure (56), 
Appendix C. 
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C. PARABOLIC GRID GENERATION 

The last grid generation technique based on the solution of P.D.E.s is the parabolic 
grid generation method. The parabolic grid generation techniques may be constructed 
by modifying elliptic methods and hence carries various advantages of the method. The 
most popular modification is the elimination of the second derivatives. The solution is 
generated by marching out in one direction like in the hyperbolic method, but the 
marching is influenced somewhat by the other boundary. Control functions can be used 
to enhance orthogonality, which would not occur normally. Because of the effect of the 
other boundaries these methods tend to have more smoothing effects than a hyperbolic 
method. The parabolic method has the characteristics that are present in both elliptic 
and hyperbolic grid generation methods. The complexity tends to be less than the el- 
liptical method and hence is faster.[Ref. 4: pp. 277-278} 
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V. RESULTS AND DISCUSSIONS 


Modern fighter aircraft designs take advantage of the strakes to improve 
controllability and enhance lift capabilities at high angles of attack. Existing aircraft are 
currently structurally modified by adding leading edge extensions to wings in order to 
improve the fuselage flow field characteristics which would eventually lead to improved 
lift and maneuverability characteristics [Ref. 28] . Complex fluid dynamic phenomena 
are associated with the high angle of attack flow over the forebody, the strake and the 
wings. The flow separates to form vortices which provide nonlinear lift. At high angles 
of attack the forebody, the strake and the wing vortices interact with each other and as 
a result self-excited unsteady flow may be triggered. When the angle of attack is further 
increased vortex breakdown occurs which will enhance flow unsteadiness and may result 
in a loss of controllability and other undesirable effects such as wing rock or tail buffet. 
Many of these interesting flow phenomena can be observed for the flow over the 
Strake-delta wing configuration model for which the grid was generated. 

In this chapter a survey of the grid generation will be done and the results of the 
numerical solution showing the characteristics of the flow field will be presented. Dis- 
cussions on vortex formation, interaction and breakdown will be made for the various 
angles of attack. 


A. GRID GENERATION 

The generation of the surface and the field grid is a prerequisite for the numerical 
solution. The grid generation can be a very time consuming process. However, the grid 
quality will determine the accuracy of the numerical solution. The double-delta wing 
analyzed here is a simple model of a modern fighter aircraft planform. The surface de- 
finition can be done entirely using linear relationships. This allowed to use relatively 
simple algebraic and geometric relationships to generate a surface grid. Even with these 
simplifications the time expended on creating a surface grid and two types of field grids 
was quite long. The amount of effort that is expended on grid generation for an actual 
aircraft configuration would very likely be more than one year. 

For symmetric bodies it is sufficient to generate half the surface grid. The gener- 
ation of the surface grid was simplified by dividing the wing into similar sections and 
writing computer programs specific for the cross-section. It was found to be simpler to 
progress from the nose to the tail by a Ax increment and compute the points (in the 














yz-plane) that represent the wing cross-section at that particular x-location. The major 


concern during this phase was to write a code that would permit a variable distribution 
of the grid lines. It was also equally important that che distribution of grid lines must 
transition as evenly as possible between fine distributions and coarser distributions. Fine 
distributions occur at locations where large numbers of grid points are required. Coarse 
distributions were chosen at locations where small grid density suffices to capture the 
physics of the flow. While a finished surface grid may appear to have smooth distrib- 
uuons and be very uniform, only the subsequent field grid generation revealed whether 
the chosen distributions were adequate. Therefore, great emphasis should be placed on 
surface grid versatility in the early stages. 

The distribution of the grid lines and the number of grid lines in the x-direction was 
changed during the interactive process of improving the surface grid to produce a suit- 
able field grid. In the case of the cylindrical grid generation, this trial and error process 
was relatively short because the field grid is composed of two-dimensional grids. On the 
other hand, for spherical grid generation both axial and circumferential surface grid dis- 
tributions were much more sensitive because they must be suitable for the generation 
of a three-dimensional mesh. During the process of generating both the cylindrical and 
spherical type grids, various small programs were written to refine and improve a par- 
ticular grid that was generated. Constructing a grid requires a knowledge and familiarity 
with grid generation codes and some prior knowledge of the flow characteristics. For 
example, the cylindrical grid required the deletion of the first points of the original sur- 
face grid in order to eliminate the singularity at the apex. Another program was then 
Written to extend the grid upstream to where conditions of the freestream where ex- 
pected. These programs can be found in Appendix F. The generation of the spherical 
grid required the writing of additional small programs in an attempt to algebraically 
modify regions of the field grid that had small discontinuities, see Appendix F. Diligence 
in changing the distributions and alteration of the nose region resulted in a smooth 
spherical field grid. It is believed that despite the larger amount of time spent for the 
generation of the spherical grid, a better quality grid compared with the cylindrical grid 
was obtained. However, because of a time constraint this could not be verified by ob- 
taining a solution on the spherical grid. The cylindrical grid was used for the numerical 
implementation because a flow solution was desired for the presentation of this work. 

A visual comparison of the two completed field grids reinforces the opinion that a 


spherical grid topology is more suitable for the subsequent numerical solution. It is 
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emphasized however, that spherical grid generation is very sensitive to the surface grid 
distributions and smoothness and requires larger computing times. 


B. FLOW FIELD CHARACTERISTICS 
The main characteristic of the flow over delta wings is the presence of the leading 
edge vortex. The nonlinear induced lift by the leading edge vortices has been actively 
investigated in recent vears. Vortical flow is an advantageous lift generation mechanism 
that can be utilized successfully at medium to high speeds. A description of the 
leeward-side flow field characteristics will be presented and the flow field structure over 
a double-delta wing at various angles of attack will be analyzed. Available experiment! 
results will be used to validate the computed results. 
1. Vortex Characteristics 

The leading edge vortices result from the roll-up of the shear layer that is shed 
from the leading edge. At moderate to high angles of attack the wingward and leeward 
side flow merges to form a shear layer that rolls up and forms a rotational vortex core. 
In the case of the double-delta wing a system of two primary vortices is formed. The 
first is along the strake and the second along the leadiriz edge of the wing section. The 
primary strake vortex reattaches itself to the centerline of the planform. The vortex 
strength is increased by a continuous feeding of vorticity from the shear lavers of the 
leading edge. Surface pressure suction peaks are produced at a location below the po- 
sition of the vortex core. Because the vortex core exhibits large gradients of vorticity 
and circumferential velocity, large viscosity effects are expected. As the vortex strength 
increases downstream. so do the lateral velocities that are near the surface. Coincidental 
with large velocities is a decrease in pressure. A secondary separation and formation of 
the secondary vortex is the result of these large lateral velocities and the associated ad- 
verse pressure gradients. If the secondary vortex is strong enough, a tertiary vortex can 
form under the secondary vortex by the same mechanisms [Ref. 17]. A schematic 
showing the leeward side vortex system and sense of rotation of these vortices is shown 
in Figure (6). Separated flow of the secondary vortex system reattaches again on the 
wing leeward surface. This separation and reattachment process can be best visualized 
with the use of surface oil flow patterns, see Figure (57). This visualization method has 
effectively shown the separation and reattachment of both primary, secondary and ter- 
tiary vortices. 

The strength of the leading edge vortex increases downstream and as the angle 
of attack is increased. The pressure gradient in the direction of a vortex core accelerate 
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Figure 6. Primary, Secondary and Tertiary Vortices 


the fluid particles until a critical angle of attack is reached. At this angle the organized 


vortex core suddenly breaks down due to the adverse pressure gradient at the trailing 
edge. . This sudden transition is more commonly known as vortex burst or vortex 
breakdown. Vortex burst is a flow phenomenon that needs to be understood because a 
loss of the suction peak will occur and this change of the induced lift can result in un- 
desirable effects on the aircraft. Work in this area by Sarpkava, Thomas, Kjelgaard and 
Sellers, Ekaterinaris, Hawk, Barnett and O’Niel, Kegelman and Roos and many others 
has been conducted in recent vears (Ref. 12,29,30,20,31,32]. 











Vortex burst is a transition from a jet-like spiraling flow to a wake-like flow. 


Adverse pressure gradient and swirl angle of the flow contribute to this phenomenon. 


Swirl angle in defined as; 
gatan ($) 


where w is the axial velocity component and y, is the azimuthial velocity component. 
Vortex breakdown will usually occur when the swirl angle exceeds a critical value of 
approximately 40°. The swirl angle and the adverse pressure gradient determine the 
type of burst for a cylindrical vortex, i.e., bubble breakdown, spiral breakdown or double 
helix breakdown [Ref. 12]. Normally, as the angle of attack is increased the burst lo- 
cation will move upstream. If the angle of attack is increased even more, wake type of 
flow behind a bluff body will be encountered. Most of the initial investigations of vortex 
generation, induced lift and vortex breakdown was completed using a single-delta wing 
configuration. In this investigation the complex flow field that results due to the pres- 
ence of multiple vortices is even further complicated for the case of the double-delta 
wing. This is due to the interaction of the strake vortex with the wing vortex and with 
the surface of the wing. The characteristics are shown for various angles of attack and 
results are discussed below in more detail. 
2. Double-Delta Wing Flow Characteristics 
a. Angle of Attack - 10° 

The computed surface flow pattern at a = 10° shown in Figure (57) does not 
indicate tertiary separation on the strake, while secondary and tertiary separation are 
shown on the wing. The leeward side flow characteristics of the double-delta wing show 
at 10.0° angle of attack are shown in Figure (58) and Figure (59). Two primary vortices 
are formed, the first is formed by the sharp leading edge of the strake and the second 
by the leading edge of the wing. Both of these vortices are continually fed by vorticity 
as they progress downstream which increases their strength. The sense of rotation for 
the primary strake vortex and primary wing vortex is the same as can be seen in the ve- 
locity vector diagrams in Figure (60) through Figure (62). Also clearly shown in these 
figures, are the primary and secondary vortices having opposite swirling directions. 
Depending on the angle of attack, both primary vortices may swirl around each other, 
but here the vortices remain separated. On the other hand, the wing tip vortex and the 


wing vortex do eventually merge as can be seen in Figure (58). 
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At 10° angle of attack the primary strake vortex separates and reattaches 


at the centerline of the wing, see Figure (57). It can also be seen that a secondary sep- 
aration occurs near the leading edge of the strake. This secondary vortex also reattaches 
itself, but it is not strong enough at 10° to generate a tertiary vortex. The wing vortex 
and reattachment of the primary and secondary vortex can also be seen in Figure (57). 
b. Angle of Attack - 19° 

At this angle of attack the primary strake vortex is much stronger. This can 
be deduced by the presence of a tertiary vortex as seen in Figure (63). This increase of 
vortex strength would produce an increase of the induced lift. The wing vortex that de- 
velops is continually fed by two sources. One source is the shear layer that is connected 
to the wing leading edge and the other is the shear layer associated with the primary 
strake vortex. This relinquishment of vorticity by the primary strake vortex causes its 
strength downstream of the kink to remain constant or even to reduce. These two 
vortices eventually merge close to the trailing edge of the wing, see Figure (64). The 
vortex burst defines the limit of vortex strength that can be maintained by the flow field. 
The burst appears to occur shortly after the primary strake and primary wing vortices 
merge, Figure (64) and Figure (65). It is at this angle of attack, i.e. just before vortex 
burst appears, where the maximum induced lift occurs. As the angle of attack is in- 
creased further the strake vortex continues to get stronger but the burst point moves 
further upstream. Close examination of particle traces, see Figure (64) and Figure (65), 
actually shows the development of a wing-tip vortex. This wing-tip will eventually merge 
with the primary wing vortex. The flow direction for both the strake and wing primary 
vortices is the same, see Figure (66) through Figure (68). These cross-sectional views 
of the velocity vectors also show the counter rotation between primary, secondary and 
tertiary vortices. At the cross-section over the wing, Figure (66), the two vortex cores 
are distinct. While at the trailing edge, Figure (68), the two vortices are merged. 

c. Angle of Attack - 22.4° 

As the angle of attack continues to increase the most prominent feature is 
probably the change in location of the vortex breakdown. The location of the burst will 
move upstream as the angle of attack is increased. Not only does the burst location 
move upstream, but the strength of the strake vortex increases. Figure (69) clearly 
shows the development of a tertiary vortex which is characteristic for a strong vortex. 
The breakdown of the vortex is readily apparent in the 22.4° flow solutions. In Figure 
(70) and Figure (71) the bursting of the strake vortex over the wing can be seen. The 
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rotations of the vortices that develop here are the same as at smaller angles of attack 
and can be verified by Figure (72) through Figure (74). 
3. Comparison with Experimental Data 

Although the purpose of Cunningham and Boer’s experiment was the investi- 
gation of unsteady phenomena, the report also presents steady state data. The devel- 
opment and location of the vortices are in qualitative agreement with Cunningham and 
Boer’s flow visualization results [Ref. 18]. These authors also present steady pressure 
measurements for five different angles of attack. Unfortunately, insufficient time was 
available toward the completion of this investigation to attempt a detailed comparison 
of the present computational results with the pressure data. The only comparison made 
was for the pressures calculated at a = 19°. Similarly, the double-delta wing studied by 
Krause and Liu [Ref. 17] and the experimental data of Brennenstuhl cited therein has 
not yet been used for comparison purposes. 

A comparison of tne pressure coefficient and spanwise location to experimental 
data can be seen in Figure (75) through Figure (77). Three axial locations on the wing 
were selected, specifically, xc = 0.40. 0.66, 0.98. The locations were selected to inves- 
tigate representative sections of the strake, the wing and the trailing edge. It can clearly 
be seen that at x.c = 0.40 there develop two suction peaks due to the primary and sec- 
ondary vortex. The location of these peaks differs from the experiment, but the trend 
is Well represented. Since the calculations are believed to be not yet fully converged, it 
is expected that more computational time will vield even more accurate results. It must 
be pointed out that these calculations are for vortical separated flow and the accurate 
capture of these characteristics is very difficult. In Figure (76), a cross-section of the 
wing shows the two suction peaks due to the primary strake and primary wing vortices 
quite well. As before. the calculations reproduce the trends quite well, but it is expected 
that more fully converged results will produce still better agreement. The last compar- 
ison is made at the trailing edge and is shown in Figure (77). Again, the trend is re- 
produced well]: but, because at this axial location vortex breakdown occurs, it is very 
difficult to achieve exact results. This graph shows a smaller pressure coefficient com- 
pared to the other axial locations. This decrease is a direct result of the vortex break- 
down. 

All the comparisons of the coefficient of pressure depicted quite well the trends 
associated with separated vortical flow. It is expected that the locations of the suction 
peaks will be even more accurate after more convergence. The comparison made here 
are based on results that required 35-40 hours of CPU time on the Crav-YMP computer. 
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It is expected that another ten hours of CPU time would produce even more accurate 
results. The expected total CPU time for each angle of attack is approximately 50 hours. 
Due to insufficient time, the fully converged results could not be presented in this paper. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


An investigation of the flow characteristics that are created in the fluid domain sur- 
rounding a double-delta wing at various angles of attack by a numerical approach was 
presented. The numerical solution of the Navier-Stokes equations requires the 
discretization of the flow domain by a smooth computational grid. Accurate represen- 
tation of the flow physics by the grid points will directly affect the quality of the flow 
field solution. Therefore, the surface and field grid density and topology must be care- 
fully chosen. With the flow field domain defined the numerical solution (in this case fi- 
nite difference) can be implemented to investigate the flow characteristics or compare 
with experimental results. 

An algebraic method of grid generation was selected for the defining of the surface 
grid and the source code is presented in Appendix E. The important precaution is that 
when developing the source code attempt to allow the flexibility of as many parameters 
as possible. The grid line distribution and number of grid points in a specified direction 
were found to be most important. Subsequent generation of the field grid will require 
the moving of grid lines to a distribution that will produce a smooth and continuous 
grid. 

Different types of numerical grid generation techniques are available, such as 
hyperbolic, elliptic and parabolic. The advantages, disadvantages and characteristics of 
these methods were discussed. The hyperbolic grid generation technique was chosen and 
two field grid topologies were generated, a cylindrical grid and a spherical grid. The cy- 
lindrical grid was easier to generate, but the spherical grid yielded a smoother grid dis- 
tribution in space. This was achieved at the expense of time and computational effort. 
However, use of the cylindrical grid would allow the generation of an acceptable flow 
field solution. 

Once the grid was completed, a finite difference algorithm in conjunction with an 
algebraic turbulence model was utilized to obtain flow field solutions at 
a = 10.0°, 19.0°, 22.4° angles of attack. Investigations of vortex generation, vortex 
interaction and vortex breakdown were conducted. At moderate angles of attack the 
double-delta wing configuration showed primary vortices generated from both the strake 
and wing. These vortices produce nonjinear vortical lift which can be very beneficial to 
fighter type aircraft that operate at high angles of attack. At approximately 19° angle 
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of attack vortex bursting occurred just after the primary strake vortex and primary wing 
vortex merged together. The formation of the primary and secondary vortices over the 
double-delta wing compared favorably with the flow visualization data of Cunningham 
and Boer [Ref. 18]. It is strongly recommended, as the next phase of this investigation, 
to compare the present numerical results with the steady pressure data obtained by these 
two authors. Using the source code presented here as a building block, future studies 
could repeat the calculations for the double-delta wing studied by Krause and Liu [Ref. 
17] and then compare with the experimental results of Brennenstuhl cited in Reference 
17. 

Future studies of this phenomenon might be to continue this same analysis utilizing 
the spherical grid to determine if the results are more accurate or if computational time 
is less, ic. the solution field converges faster. Furthermore with the recent increased 
interest in dynamic stall phenomenon, an analysis could be done to compare the com- 
putational results of a pitching straked-wing to the experimental studies done by 
Cunningham and Boer. Because a major portion of time is spent in the generation of a 
field grid for an analysis of this type, the work load would be reduced by using the grid 
presented in this thesis. 

The work load involved in generating a field grid is significantly increased when the 
body under investigation is an entire aircraft. For this reason the need for a method of 
quickly producing a surface grid would expedite a numerically generated solution of a 
flow field and allow more time for the improvement of numerical methods. 
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APPENDIX A. SURFACE GRID FIGURES 
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Double-Delta Wing Surface Grid (120x240x 1) 


Figure 8. 
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Double-Delta Wing - top view (100x240x1) 


Figure 9. 
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Figure 18. Grid Distribution of the Wing (26x240x 1!) 









































Figure 18. Grid Distribution of the Rectangular Section (16x240x1) 
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Figure 23. Cylindrical (H-O) Grid Topology (130x240x68) 
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Figure 24. Spherical (C-O) Grid Topology (160X240x68) 


70 








APPENDIX B. FIELD GRID FIGURES -- CYLINDRICAL 
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Cylindrical Grid Configuration (130x240x68) 


Figure 25, 
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Figure 27. Cylindrical Grid - side view (130x240x68) __ 
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Figure 36. Typical Cross-section of the Rectangular Section - front view 
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Figure 39. Typical Cross-section of the Wake - front view 
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APPENDIX C. FIELD GRID FIGURES -- SPHERICAL 
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Spherical Grid Configuration (160x240x68) 





Figure 42. 
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Spherical Grid - side view (160x240x68) 


Figure 44. 














PAU 


i\\\\ 
oor 


93 


HTT 
MYL 
HD mL 







WW 










Typical Cross-section of the strake - front view 


Figure 45. 
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Figure 48. Typical Cross-section of the Wing - front view 
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Figure $1. Typical Cross-section of the Rectangular Section - front view 
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Figure 54. Typical Cross-section of the Wake - front view 
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APPENDIX D. RESULTS AND DISCUSSION FIGURES 








105 














(gOXE9XDL) ‘SAOIXS'E =F “7Z7'O=W - oOl 3B W4azeq MOLY ad"JINS = *7g aundry 





106 








(S9XE9XOL) “OAOTX8'E =2e ‘770 = - oOL 38 Sadvdy ajeg = “gg andi d 














(S9XE9XOL) “OAOIXSE = 2H '7CO=W - Ol 28 UONEIOT Xap0,A 





"6S andy 


108 








9 andy 


\ 


QAOTX8'E = 2a ‘7770 =W - oOT JE S40999A AYID0J9Q aYENS — *g 


= 
4 


Mtl 


nit 


Se a, 


a il OG pe 


Te 


aii tt 






nih! 


PT 





ee 





OAL 
THLE Ahy 


a 


\ < era 4//// 
ye MS 1/1) 





\ 





| 

\ 
\ 
BS 


\ 


\ 
are 
\ 


; 
; 
| 


TULA ys! 


109 














\ 


\ 


2 ee Oe SS 
_ _ -_ a ee 
_ a ee — = 
“ ager -- 
was A “Loca 
cS we 4 4A ALE 
gE. fe ee hs 
a a LILI 
ye 
ae ese 
- v - 2 / / | 
Y “ Py; 7 / / i | \ 
Bo ee Fh TN ANNS 
fe 90 opp FANART 
/ 
5 A FELIU 
/ - epee) 
7 a EY at 
ye fy EDP TOL 
R\ 
, / ‘ / 44s fheer i 
atti 
. / / / / ft? ff fee 
7 / ‘ / | | Te oa ea as 
, / /- { { \ VAS 
F ‘ ‘ \ \ \ V\\NH? 


NNN? 


110 

























MIINVAN SoS 
HHUA ANN AQ ~~ 


BAKKE SNS 


AWA SN NNN 


AMIS 
MWA SS NN 


‘\ 


“ 


WWWANNN NNN NON OY \ 


EEE Se 2 XN 


MUAY SR QO Y 


BUAVAN ANN N ON NNN XN 


MAAWAL ARNON ON NON OON ~ 


SSS SSS Nn ~ 


MAAR LAN NS NS NON ON ~ 


MAAN SS NN ONS ~ > 


“N 
MINNA QQ 
‘ “N 
MMA NN QQ QL 
“ 
BIWWVVN NNN NY QS 


3.8x LNE6 


Wing Velocity Vectors at 10° - M/=0.22, Re 


Figure 61. 

















ee 

















aS 
~ 
\ al 
BWA NNN 
AW 
PBI 
BOY LNA NNN NON 
\ 
MWA SN SSN 
\ Sm 
PNA ae SR NRE 
\ “NX 
WWW Ge Phage ae 
. : Be 
PE NSN Ard LO eS 
" ~ ie 
BMOMUV A 7 ge 






| MWA \ \ Nv KY 






HE PBA NON 





3.8x 1OE6 


 caaenaciaaa NNN 


















BAAN NSN NN ON ONO x 














SSSANS SNS SS NON ~ ~ ~ ~ 





Ah We aoe 






MA fl 


‘ ' t \ \ VV NNN | 


LU A 


7 
‘ t \ \ \ VON ANSE 


T. E. Velocity Vectors at 10° - M[= 0.22, Re 





Crea! Wea: ee Year SS) p 2 ~ ~ 









eee 7/7 


Figure 62. 
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Figure 63. 
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Figure 67. 








COONS ie 


es 


AWWUIAN NL 


ag! e 
AMMAN NNN \ 


ae WS AS 


MWA NAN NN 


3 SSS OO 


| 
Sur Rs ll | 
foes —_ HMM 
So HW 


So US esl 7 


117 


H . 
| oe NNNONON ~ 


saNANN NNN ~ ~ ~ ~ 


SR i ia = “— 





3.8x 10E6 


0.22, Re= 


T. E. Velocity Vectors at 19° -M 


Figure 68. 








Sif 
=| (i 
——=—= 


(SOXE9NOL) “9AOIX8'E = 2Y ‘770 =I - ob TT JE Udayed MOPy adeying = “6g aang q 





118 








(S9XE9XOL) “9AOIX8'E = 94 “770 =A - ob ZZ FE Savy ape «OL aanaiy 


119 











(S9XE9XOL) ‘YFOIX8'E =F “ZT = IN - of TZ JE UONLIOT XBHOA 





"TL aaingiy 


120 














9301X8" 
€=9u° 
770 = 
= W - 
ard . 
TZ JE S10j93A_ A}IDOI9 
ID0J9A, 9YvI 
~  * 
TL aanai 
7 14 


= 


mae 





till 


Ma 


S= 
SS 
=== 





SS | 
=—— = 
SS 
~ ces 
SS 
eS 
Se © 2 
ye est 
\ 
Vee 
Mya g 
ph 
bp yo 
| / 
a aa, | 
f/f 








—_ 
a 
— 























BWW 
BIW A SD 
AWWA A 
IWAN NN NN 


3.8x 10E6 





0.22, Re 


NAAN ANN NNO 





-M 


MIWA NNN NNO 
AWWA ANNOY 
AWWA NNN 


BIA N NNN NONOONON 


MWA NNN NNN 


Wing Velocity Vectors at 22.4° 





SON Nee eld EOLA KR NNO 






~ ~ 
7 / i aN \ NON SR ee ee MAMAN AN SRN NS 


Figure 73. 




















eos" 
yw 


YVAN 
VV ANN NAA 
VV AUAAAAN 
VV AVA 
4 oe VV AVAAN 


Ve VAUAALAN 


PoP A © A NAAAAAN 


123 





\ 


=. 
as 
RRQLOOS 
WSS 
SQQancdk Ss 
{ Wan eat 
A Qa NS 
NWWQA QS \ 
BNW ON SS 
BWW OND 


oS 
SN 
SS 





™~. 

~ 
BAIA A = 
™~ 


AIAN AA ae 


BINNS SN 
™~ 


BION NNN NNN ~N 


OSS a Ge ~N 


MMA ARR ~ 


EIA ~™ 


i i i ~ ~ 


ee ee ~~ _ 





3.8x 10E6 


0.22, Re= 


T. E. Velocity Vectors at 22.4° - M 


Figure 74, 














Ob’ = 3/X 3e JUAaIIYJa0-) auNnsSadg adeNG “SL aan3ly 


q/A uones07 asimueds 
8'0 90 v0 | 


Or'0=2/X pamseayy @ 
Or'0 = 9/x pandumy ——}-+t 


"O 
Lm 4 
4") 
” 
” 
ce 
a 
te) 
C) 
oO 
@ 
=a 
= 
2H 
© 
=| 
cr 
0 
C-) 
x) 


v 


gO X SE = 24.61 = ?'TO=W 
Ov'O = 2/x ie JUTIDYJIO7) JANSSIIdg BIeEpNS 








99°Q = 3/X JU JU JIOD ainssaig adeying gf andy 


q/A uoije907 asimueds 
80 90 r'0 re 


125 


99°0 = 9/x pamsvoyy i 





"O 
Lm 4 
o 
~ 
” 
c 
= 
ro) 
‘2 
° 
rc) 
= 
2. 
© 
S 
er 
99°0 = 9/x pandwoy —— ; 
= : C) 
Tc 


=z 


90L X SE=29Y‘ 6T= 72 TO=W 
99°Q = 9/x 1e JUAIDYJa07 aANSsaig adeyNS 











86) = 3/X ye JUsId aor ainssaig advpUNG *Z/ sandy 





q/A uones07 asimueds 
01 80 9°0 v0 70 00 






© 


— 


N 


860 =5/x pamsssyy A | 
860 = 9/x panduwoy —— 


foe) 





dy - ywoanyya0y aunssaig 


=z 


90T X SE = 2Y 61 =”'TO=W 
96'0 = 2/X Je JWwUAIDIja0D aanssaig aIeyING 














APPENDIX E. SOURCE CODE FOR SURFACE GRIDS 


127 























NAANARNANANANAANANANANAANNANAANNNNAN aNNNNANANNANANNNANNANANNNNNANNNNANANNANANNA 


THIS IS A PROGRAM TO GENERATE A SURFACE GRID FOR A DOUBLE DELTA 
WING AIRFOIL THAT WILL BE ANALYZED FOR MY MASTERS THESIS. 


CERTAIN DATA WILL BE REQUIRED VIA A DATA FILE AND INCLUDES THE 
FOLLOWING VARIABLES: 


NOSERAD - THE RADIUS DECSCRIBING THE NOSE 
NOSGDX1 - THE NUMBER OF GRID POINTS, IN THE X DIRECTION DESIRED 
IN THE ROUNDED NOSE 


NOSGDX2 - THE NUMBER OF GRID POINTS, IN THE X DIRECTION DESIRED 
IN THE FIRST DELTA WING. DO NOT INCLUDE LAST GRID IN 
THE NOSE ROUNDING 

NOSGDY - THE NUMBER OF GRID POINTS, IN THE Y DIRECTION DESIRED 


IN THE NOSE BODY 
CRVGDY - THE NUMBER OF GRID POINTS, IN THE Y DIRECTION DESIRED 
IN THE ROUNDED LEADING EDGE 


LEN1 - THE LENGTH FROM THE LEADING EDGE TO THE SECOND WING 
LEN2 - THE HALF WIDTH OF THE FIRST DELTA WING 
LEN3 - THE HALF THICKNESS OF THE FIRST DELTA WING 


THE FOLLOWING IS A DEFINITION OF SOME OF THE VARIABLES USED IN 
THE FIRST PART OF THE PROGRAM: 


DELANG THE DELTA ANGLE USED TO GENERATE GRIDS IN THE Y DIRECTION 
DISTA - THE DISTANCE FROM THE TIP OF THE AIRFOIL TO THE LAST 
GRID GENERATED IN THE NOSE ROUNDING 


DELDIST - THE DELTA DISTANCE IN THE X DIRECTION USED IN ROUNDING 
THE NOSE 

ANG - THE ANGLE USED TO CALCULATE THE SPECIFIC GRIDS IN THE 
Y DIRECTION 

RAD - THE RADIAL DISTANCE USED TO CALCULATE GRIDS IN THE Y 
DIRECTION 

Xl - THE DISTANCE FROM THE CENTER OF THE NOSE TO THE TRAILING 


EDGE OF THE FIRST DELTA WING 

ANG1 - THE ANGLE BETWEEN THE CENTERLINE AND TRAILING EDGE 

ANG2 - THE ANGLE BETWEEN THE TRAILING EDGE AND LEADING EDGE 
OF THE FIRST DELTA WING 

THETAl - THE ANGLE FORMED BY THE NOSE ROUNDING, IT WILL BE 
PERPENDICULAR TO THE LEADING EDGE OF THE FIRST DELTA WING 


DISTB - THE DISTANCE FROM THE TRAILING EDGE TO THE NOSE RADIUS 

INUM - THE TOTAL NUMBER OF LINEAR GRIDS ON THE UPPER AND LOWER 
SURFACE 

DISTF - THE DISTANCE IN THE 2 DIRECTION TRAVERSED BY THE THICKNESS 


OF THE LEADING EDGE 
THETA2 - THE ANGLE FORMED BY THE THICKNESS OF THE FIRST DELTA WING 


DIStTc - THE HALF WIDTH IN THE Y DIRECTION USED TO GENERATE Y AND 
Z DIRECTION GRIDS 

DISTD - THE HALF THICKNESS USED IN SAME CALCULATIONS AS DISTC 

THETA3 - SAME TYPE OF ANGLE AS THETAl BUT USED FOR Y AND Z GRIDS 

DISTE - THE DISTANCE IN THE Y DIRECTION FOR Y AND Z GRID 

DELY - THE DELTA Y DISTANCE FOR THE DISTE VARIABLE 


DIMENSION X(170,150,1),¥(170,150,1) ,2(170,150,1) 

INTEGER NOSGDX1, NOSGDX2 , NOSGDY , CRVGDY , BODGDX1, BODGDX2 , NOSGDX, 
7 SIDGRD, WAKGPD, CRVGDY1 : 
REAL LEN], LEN2, LEN3,LEN4, LENS, LEN6, NOSERAD, NOSERAD1 











OPEN (UNIT=10, FILE='ddwsg.in',STATUS='OLD') 
c OPEN (UNIT=12, FILE='ddwsg.dat',STATUS='NEW') 


READ(10,*) NOSERAD, NOSGDX1,NOSGDX4 , NOSGDX5 , NOSGDY 
READ(10,*) CRVGDY,SIDGRD,WAKGRD 

READ(10,*) BODGDX1,BODGDX2,THETA4 

READ(10,*) LEN1,LEN2,LEN3 

READ(10,*) LEN4, LENS, LEN6 

READ(10,*) DISTX1A,DISTX1B, DISTX2, DISTX3 , DISTX4 
READ(10,*) C1A,E1A,C1B,E1B,C2A,E2A,C2B,E2B 
READ(10,*) C3A,E3A,C3B,E3B 

READ(10,*) WAKLEN 


PRINT *,*YES1'* 


c 
C THIS SECTION DOES SOME PRELIMINARY CALCULATIONS ON THE WING 
Cc 
PI=2.0*ASIN(1.0) 
DELANG=P1/ (NOSGDY-1) 
X1=( (LEN1=NOSERAD) **2.0+LEN2*#*2.0) #*0.5 
DUMMY=LEN2/X1 
ANG1=ASIN (DUMMY) 
DUMMY=NOSERAD/X1 
ANG2=ACOS (DUMMY) 
3 THETA1=PI-ANG1-ANG2 
c 
C THIS SECTICKX DOES THE ROUNDING OF THE NOSE TO AVOID SINGULARITIES 
c 


DO 10 IX=1,NOSGDX1 
DO 10 IY=1,NOSGDY 
IF (IX.EQ.1) THEN 
X(IX,IY¥,1)=0.0 
¥(IX,IY,1)=0.0 
2(IX,I¥,1)=0.0 
ELSE 
X (IX, IY, 1) =NOSERAD-COS ( (THETA1/ (NOSGDX1-1) ) * (IX=2) ) 
* *NOSERAD . 
ANG=PI/2.0+(IY-1) *DELANG 
RAD=NOSERAD*SIN ( (THETA1/ (NOSGDX1-1) ) *(IX-1) ) 
¥(IX,IY,1)=-1.0*COS (ANG) *RAD 
2(IX,IY,1)=SIN (ANG) *RAD 
ENDIF 
10 CONTINUE 


PRINT *, 'YES2' 


c 
C THIS SECTION DOES SOME PRELIMINARY CALCULATIONS ON THE WING 
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DISTB=LEN1~NOSERAD+COS (THETA1) *NOSERAD 
INUM=NOSGDY-CRVGDY 
DISTF=LEN3~-SIN (THETA1} #NOSERAD 
THETA2=ATAN (DISTF/DISTA; 

CUMDIST=0.0 


THIS SECTION COMPUTES GRID POINTS FOR THE NOSE UP TO THE SECOND WING 


CRVGDY1=NOSGDY-104 
NOSGDX2=NOSGDX4/2+NOSGDX5/2 


DO 20 IX=1,NOSGDX2 


IF(IX.LT.4 ) CRVGDYl=CRVGDY1-2 
IF(IX.GE.4.AND.IX.LT.6) CRVGDY1=CRVGDYi-2 
IF(IX.GE.6.AND.IX.LT.NOSGDX2) CRVGDY1=CRVGDY1-2 
IF (CRVGDY1.LT.7) CRVGDY1=7 

INUMN=NOSGDY-CRVGDY1 

CUMDELY=0.0 


IF (IX. LE.NOSGDX4/2) THEN 

CALL STRCH4 (DISTB, NOSGDX4 , IX, DELDIST, DISTX1A) 
ELSE 

CALL STRCH4 (DISTB, NOSGDX5 , 1X, DELDIST , DISTX18) 
ENDIF 


ICOUNT=0 
CUMDIST=CUMDIST+DELDIST 

DISTC=SIN (THETA1) #*NOSERAD+ (CUMDIST) /TAN (THETA) 
DISTD=SIN (THETA1) *NOSERAD+ (CUMDIST) *TAN (THETA2) 
X1= ( (DISTC-NOSERAD) **2.0+DISTD##2.0) **0.5 
DUMMY=DISTD/X1 

ANG1=ASIN (DUMMY) 

DUMMY=NOSERAD/X1 

ANG2=ACOS (DUMMY) 

THETA3=PI-ANG1-ANG2 
DISTE=DISTC-NOSERAD+COS (THETA3) *NOSERAD 
DELANG=2. 0*THETA3/ (CRVGDY1+1) 


DO 20 IY<1,NoSGDy 
X (IX+NOSGDX1, IY, 1) =X(IX+NOSGDX1-1, IY, 1)+DELDIST 
IF (IY. LE. (INUMN/2) ) THEN 


¥ (IX+NOSGDX1, IY, 1) =CUMDELY 
2 (IX+NOSGDX1, IY, 1) =DISTD- (CUMDELY/TAN (THETA3) ) 


CALL STRCH9 (DISTE, INUMN/2-1,IY,DELY,C1A,C1B,EJA,E1B, 
* NOSGDX2 , IX) 


CUMDELY=CUMDELY+DELY 
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ELSEIF (IY.GE. (INUMN/2)+1.AND.IY.LE. (INUMN/2) +CRVGDY1+1) THEN 
ICOUNT=ICOUNT+1 
DUMMY=THETA3-ICOUNT*DELANG 
¥ (IX+NOSGDX1, TY, 1) = (DISTC-NOSERAD) +COS (DUMMY) *NOSERAD 
2 (IX+NOSGD%1,1¥,1)=SIN (DUMMY) *NOSERAD 
ELSE 


Y (IX+NOSGDX1, IY, 1) =¥ (IX+NOSGDX1 , NOSGDY+1-IY,1) 
Z (IX+NOSGDX1,IY,1)=-2 (IX+NOSGDX1, NOSGDY+1-IY,1) 


ENDIF 
CONTINUE 


PRINT *, 'YES3! 


BELOW ARE LISTED SOME MORE VARIABLES USED IN THE PROGRAM: 


BODGDX1 + THE NUMBER OF SECOND DELTA WING GRIDS UP TO THE RECTANGULAR 
SECTION AND NOT INCLUDING THE FIRST GRID FROM THE NOSE 

NOSGDX -~- THE TOTAL NUMBER OF GRIDS IN THE X DIRECTION OF THE FIRST 
DELTA WING 

THETA¢ - THE ANGLE FORMED BY THE SECOND DELTA WING 


DISTZU1 - THE DISTANCE IN THE POSITIVE Z DIRECTION FOR THE FIRST NACA 
CROSS-SECTION ENCOUNTERED (DISTZU) 

DISTZD1 - THE SAME DISTANCE ON THE LOWER SURFACE (DISTZD) 

LEN4 - THE LENGTH IN THE X DIRECTION FROM THE SECOND DELTA WING TO 
THE REAR RECTANGULAR SECTION 

LENS - THE TOTAL LENGTH OF THE SECOND DELTA WING 


THETAS = THE ANGLE ON THE UPPER SURFACE FROM THE CENTERLINE TO THE 
FIRST NACA SECTION ENCOUNTERED 


THETA6 ~ THE SAME ANGLE ON THE LOWER SURFACE 

DISTA - THE DISTANCE IN THE Y DIRECTION WITH A NACA CROSS-SECTION 
CHRD ~ CHORD LENGTH OF THE NACA SECTION 

XDIST ~ THE POSITION IN THE X DIRECTION ON THE NACA AIRFOIL 


THE SECTION BELOW GENERATES THE GRID FOR THE SECOND DELTA WING 
NOSGDX=NOSGDX1+NOSGDX2 
CUMDIST=0.0 
DO 30 IX=NOSGDX+1 , NOSGDX+BODGDX1 
IXI=1X-NOSGDX1-NOSGDX2 
CALL STRCH4 (LEN4 , BODGDX1, 1X1, DELDIST, DISTX2) 
CUMDIST=CUMDIST+DELDIST 
X(IX,1,1)=X(IX-1,1,1)+DELDIST 
Y¥(IX,1,1)=¥(IX-1,1,1) 
Z(IX,1,1)#Z(1X-1,1,1) 


DISTE=LEN2+ (CUMDIST/TAN (THETA4) ) “NOSERAD/SIN (THETA4) 
ICOUNT=0 
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JCOUNT=0 
KCOUNT=0 
CUMDELY=C .0 

DO 30 Iy=2,NOSGDY 


CALL STRCH8 (DISTE, INUM/2-1, I1Y-1,DELY,C2A,C2B,E2A,E2B, 
BODGDX1, IX-NOSGDX) 


CUMDELY=CUMDELY+DELY 
X(IX,I¥,1)=X(IX-1,1Y,1)+DELDIST 


IF(IY.LE. (INUM/2) ) THEN 
¥ (IX, IY,1)=¥ (IX, Iy-1,1)+DELY 


ELSEIF (IY.GE. (INUM/2)+1.AND.IY.LE. (INUM/2) + 
(CRVGDY-13) /2+1) THEN 


THETA9=ATAN ( (¥ (IX, INUM/2,1) -¥ (IX, INUM/2-1,1))/ 
(2 (IX, INUM/2-1,1)-% (IX, INUM/2,1))) 
DELANG=2.0*THETAS/ (CRVGDY+1) 
JCOUNT=JCOUNT+1 
YRAD=Z (IX, INUM/2,1) /SIN(THETA9) 
¥(IX,IY,1)=¥ (IX, INUM/2,1) + (COS (THETA9-DELANG*JCOUNT) 
*YRAD) ~YRAD*COS (THETAQ) 
ELSE 
Y¥(IX,IY,1)=¥ (IX, NOSGDY+1-IY,1) 
ENDIF 
CALL NACA006 (DISTZU1, DISTZD1, LENS, CUMDIST) 


THETAS=ATAN ( (2 (IX,1,1)-DISTZU1) /LEN2) 
THETA6=ATAN ( (Z(IX,1,1)+DISTZD1) /LEN2) 


IF (IY.LT. ( (NOSGDY-1) /2) +1) THEN 

IF (¥(IX,IY,1).LE.LEN2) THEN 
Z(IX,1Y,1)=Z(IX,1,1)-TAN(THETAS) *¥ (IX, IY,1) 

ELSEIF (¥(IX,IY,1).GT.LEN2.AND.IY¥.LE. (INUM/2) ) THEN 
DISTA=Y (IX, IY,1) -LEN2 
CHRD=LENS-TAN (THETA4) *DISTA 
XDIST=CUMDIST~TAN (THETA4) ®DISTA 
CALL NACA006 (DISTZU, DISTZD, CHRD, XDIST) 
2 (IX, I¥,1)=DISTzu 

ELSE 


KCOUNT=KCOUNT+1 ; 
Z2(IX,IY,1) =SIN (THETA9-DELANG*KCOUNT) *YRAD 


ENDIF 














ELSEIF (IY.EQ. ( (NOSGDY-1) /2)+1) THEN 


KCOUNT=KCOUNT+1 
Z2(IX,1Y,1)=0.0 


ELSE 
IF(1Y.LE.INUM/2+CRVGDY+1) THEN 


KCOUNT=KCOUNT+1 
Z(1IX, 1Y,1)=SIN (THETA9-DELANG* KCOUNT) *YRAD 


ELSEIF (¥(1X,1Y,1).GT.LEN2) THEN 
DISTA=Y (IX, IY, 1) ~LEN2 
CHRD=LEN5~-TAN (THETA4 ) *DISTA 
XDIST=CUMDIST-TAN (THETA4) *DISTA 
CALL NACA006 (DISTZU, DISTZD, CHRD, XDIST) 
Z2(IX,IY,1)=DISTZD 

ELSE 


ICOUNT=ICOUNT+1 
Z(IX,IY,1)=Z(IX-1,NOSGDY,1)+TAN (THETA6) *Y (IX, IY,1) 


ENDIF 
ENDIF 
30 CONTINUE 
PRINT *, 'YES4! 
c 
C THIS SECTION COMPUTES THE LAST RECTANGULAR SECTION OF THE WING 
Cc 
CUMDIST=0.0 
INUM2=NOSGDY-SIDGRD 
NOSERAD1=NOSERAD*2.0 
DO 40 IX=NOSGDX+BODGDX1+1 , NOSGDX+BODGDX1+BODGDX2 


IF (IX.GE.NOSGDX+BODGDX1+10) INUM2=INUM2+3 
XT=( (LENS-LEN4-NOSERAD) **#2.0+ 


* (Z (NOSGDX+BODGDX1, IY, 1) -NOSERAD) **2.0) *#0.5 
ANG1=ATAN ( (Z (NOSGDX+BODGDX1,IY,1) -NOSERAD) / 

* (LENS-LEN4-NOSERAD) ) 
ANG2=ACOS (NOSERAD/X1) 


THETA7=PI-ANG1-ANG2 

DELANG=THETA7/2.0 

DISTD= ( LEN5- LEN4 -NOSERAD+COS (THETA?) *NOSERAD) 
CUMDELY=0.0 

ICOUNT=0 

JCOUNT=0 

KCOUNT=0 

LCOUNT=0 
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MCOUNT=0 
IXI=IX-NOSGDX1-NOSGDX2-BODGDX1 


CALL STRCH4 (DISTD, BODGDX2 , 1X1, DELDIST, DISTX3) r 
CUMDIST=CUMDIST+DELDIST 
DO 40 IY=1,NOSGDY 

XI=( (LENS-LEN4-NOSERAD) **2.0+ 


* (2 (NOSGDX+BODGDX1, IY, 1) -NOSERAD) **#2.0) #*0.5 
ANG1&ATAN ( (Z (NOSGDX+BODGDX1, IY, 1) -NOSERAD) / 

* (LEN5~LEN4-NOSEFAD) ) 
ANG2=ACOS (NOSERAD/X1) 


THETA7=PI-ANG1-ANG2 
DELANG=THETA7/2.0 
DISTD=(LENS5-LEN4 -NOSERAD+COS (THETA?) *NOSERAD) 
X(IX,IY,1)=X(1IX-1,1Y,1)+DELDIST 
DISTE=LEN6-NOSERAD 
IF (IY.LE. ( (NOSGDY-1) /2)+1) THEN 

IF (IY.LE. (INUM2/2)) TREN 

¥ (IX, IY, 1) =CUMDELY 


CALL STRCH8 (DISTE, INUM2/2-1,IY,DELY,C3A,C3B,E3A,E3B, 
7 BODGDX2 , IX-NOSGDX-BODGDX1) 


CUMDELY=CUMDELY+DELY 


ELSEIF (IY.GT.INUM2/2.AND. 
* IY. LE. INUM2/2+2) THEN 


ICOURT=ICOUNT+1 
DELANG1=(PI/4.vU) 
¥ (IX, IY,1)=¥ (IX, INUM2/2,1)+ 
* COS (PI/2. 0-ICOUNT*DELANG1) *NOSERAD 

ELSE 
¥ (IX, IY,1)=¥(IX,1¥-1,1) 

ENDIF 

ELSE 


JCOUNT=JCOUNT+1 
Y¥(IX, 1¥,1)=¥(IX, 1Y-2*JCOUNT, 1) 


ENDIF 
IF(Y(IX,IY,1).LE.LEN2.AND.IY.LE.INUM2/2) THEN 
Z2(IX,IY,1)=Z(IX-1,1Y, 1) -DELDIST/TAN (THETA?) 
IF(Z(1X,1IY,1) .LT.NOSERAD1) THEN 
2(IX,I¥,1)*NOSERAD1 











ENDIF 

ELSEIF (¥(IX,1IY,1).GT.LEN2.AND.IY.LE.INUM2/2) THEN 
DISTA=Y¥ (IX, IY¥,1) ~LEN2 
CHRD=LEN5-TAN (THETA4 ) *DISTA 
XDIST=CHRD- (LEN5-LEN4) +CUMDIST 
CALL NACAO006 (DISTZU, DISTZD, CHRD, XDIST) 


Z(IX,1TY,1)=DIST2ZU 
IF(Z(IX,1IY,1).LT.NOSERAD1) THEN 


2 (IX, IY,1)=NOSERAD1 
ENDIF 
ELSEIF (IY.GT. INUM2/2.AND.1Y.LE. INUM2/2+2) THEN 
LCOUNT=LCOUNT+1 
2 (IX, 1Y,1)=2 (IX, INUM2/2, 1) -NOSERAD+ 
* SIN (PI/2-LCOUNT*DELANG1) *NOSERAD 
ELSEIF (IY.GT.INUM2/242.AND.IY.LT. (NOSGDY-1) /2+1) THEN 
DIFF=((NOSGDY-1) /2+1) - (INUM2/2+2) 
DEL=Z (IX, INUM2/2+2,1)/DIFF 
2(IX,IY,1)=Z(1IX,IY-1,1) <DEL 
ELSEIF (IY.EQ. ((NOSGDY-1) /2) +1) THEN 
2(1X,1Y,1)=0.0 
ELSE 


MCOUNT=MCOUNT+1 
Z(IX,IY,1)=-1.0%*2 (IX, 1Y-2*MCOUNT, 1) 


ENDIF 
40 CONTINUE 
PRINT *,*YESS! 


THIS SECTION DOES THE REPEATING OF THE TRAILING EDGE GRID TO 
FORM THE SECTION OF THE SURFACE GRID THAT FORMS THE WAKE 


aaan 


TOT=NOSGDX+BODGDX1+ BODGDX2 
DO 50 IX=1,WAKGRD 
CALL STRCHS5 (WAKLEN, WAKGRD, 1X, DELX, DISTX4) 
DO 50 IY=1,NOSGDY 
ITOT=TOT+IX 


X(ITOT, IY,1)=X(ITOT-1,1Y,1)+DELX 
Y¥(ITOT, IY,1)=¥(ITOT-ix,IY,1) 








Z(ITOT,IY,1)=Z2(ITOT-ix,1yY,1) 
50 CONTINUE 


PRINT *,‘YES6' 


Cc 
C PRINT THE DATA TO A FILE 
Cc 
Cc DO 11 I=40,42 
c DO 11 I=1,NOSGDX+BODGDX1+BODGDX2+WAKEGRD 
c DO 11 J=1,NOSGDY 
c PRINT *,1,3,X(I,J,1),¥(I,J,1),2(I,3,1) 
cll CONTINUE 

IZ=1 

IX=1 

LY=163 

REWIKE 3 

WRITE (3) IY-IX+1,NOSGDY,1Z 

WRITE(3) ((X(1,3,1),I=IX,I¥) ,J=1,NOSGDY) , 

* ((¥(I,3,1),I=IX, IY) ,J=1,NOSGDY) , 

* ((Z2(1,d,1),1=IX, IY) ,J=1,NOSGDY) 

CLOSE (UNIT=10) 

CLOSE (UNIT=11) 

STOP 

END 
eccceeceececceceeceeceececeeecececececececececeeecceccecceceeceeceecececeeeecec , 
Cece eeeeeeeececrec cecer cercecermececececcecerecrecarececeecenereeceececececcr 
Cc 


2 THIS IS THE SUBROUTINE NACAOO6 AND IT EVALUTATES THE Z DIRECTION VALUES c 
THAT ARE ASSIGNED TO THAT SPECIFIC CROSS-SECTION 

: C 

eceecececcceccececceccccecceeeeeceecceceececccecceececcececcececccececcceccecccce 

eccecceecececeececcececceceececcecececeeceececececceecceccececcecccceecccccccce 


SUBROUTINE NACA006 (ZPOS, ZNEG, CHORD, XVAL) 


DIMENSION XPC(26) ,YPC(26) 


Rae XPC/0.0,0.5,0.75,1.25,2.5,5.0,7.5,10.0,15.0,20.0,25.0,30.0, 
35. 0, 40.0,45.0,50.0, 55. 0,60.0,65.0,70.0,75.0,80.0,85.0, 


* 90.0,95.0,100. 0/ 

DATA YPC/0.0,0.494,0.596,0.754,1.024,1.405,1.692,1.928,2.298, 

* 2.572,2.772,2.907,2.981,2.995,2.919,2.775,2.575,2.331, 
* 2.050,1.740,1.412,1.072,0.737,0.423,0.157,0.0/ 


X=XVAL/CHORD*100.0 
DO 10 I=1,25 
IF(X.GT.XPC(I) .AND.X.LT.XPC(I+1) ) THEN : 


IVAL=I 








GO TO 99 
ELSE 
CONTINUE 
ENDIF 
10 CONTINUE 
99 VAL1=X-XPC (IVAL) 
VAL2=XPC (IVAL+1) ~XPC (IVAL) 
VAL3=YPC (IVAL+1) ~YPC(IVAL) 
ZVAL= (VAL1*VAL3/VAL2) +Y¥PC (IVAL) 
2POS=ZVAL*CHORD/100.0 
ZNEG=-1.0*ZPOS 


RETURN 


ececececececcecceeccecececcececcecececececceccecceccceccececcc 
cecceceececeececcececceccececcceccececeececececeeccececececc 


c c 
C THIS SUBROUTINE STRETCHES THE GRID IN THE X DIRECTION c 
c c 


ecceccececeeceecececececcececcecceececcececceeececececcecececcce 
ececeeceeccec\ececececceccecececcecececeeccecececececcecceccccce 


SUBROUTINE STRCHi(DIA,IDIV,IA,DELDIST) 
YDIST=DIA/IDIV 
IF(IA.LE.IDIV/2) THEN 


IB=IA 
Ic=1 


ELSE 


IB=IDIV-IA 
Ic=-1 


ENDIF 


_IF(IA.GE.IDIV/2.AND.IA.LE.IDIV/2+1) THEN 


DELDIST=((DIA/2.0) **2.0-((IDIV/2=1) #YDIST) ##2.0)#*0.5 


ELSE 
DIST1=( (DIA/2.0) **2.0-(IB*YDIST) **2.0)**0.5 
DIST2=( (DIA/2.0) **2.0-((IB-IC) *¥YDIST) **2.0) **0.5 
DELDIST=ABS (DIST1-DIST2) 

ENDIF 


RETURN 
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END 


cececcececcecceccececcceceeceeececcececceecececcceccccecccecccccce 
eccecececceccececececccececeececccececeecceccececccececececececcecc 


c Cc 
C THIS IS THE SUBROUTINE FOR STRETCHING IN THE Y¥ DIRECTION Cc 
ce c 


cecececeecececeececececececececececececeeccececceececeeceeceecec 
cececeececececeeceececeeeceeceececeeccecceeceeceececeeeceecececcce 
SUBROUTINE STRCH2(RAD,JDIV, IA, DELY) 
YDIST=RAD/JDIV 
IF (IA.EQ.1) THES 
DELY=((RAD**2.0)-((JDIV-IA) *YDIST) #*2.0) **0.5 
ELSE 
DIST1=((RAD**2.0)-((JDIV-IA) *YDIST) **2.0) *#0.5 
DIST2=((RAD**2.0)-((JDIV-IA+1) *¥YDIST) **2.0) ##0.5 
DELY=ABS (DIST1-DIST2) 
ENDIF 
RETURN 
END 
ecceceeeececeecececececececceccececececcececececeececececcececec 
ececeeceeececcecceceeceececececeeecceceececececceececececececccc 
Cc Cc 
C THIS SUBROCTINE 1S FOR LINEAR STRETCHING IN THE X DIRECTION C 
c Cc 
ececeeecceccececeeceececeeceeeceecececeececcececececececececcecce 
ceeeceeeececcceecceceececceceececeeceececceececeececececeecceccece 
SUBROUTINE STRCH4 (DIA, IDIV, IA, DELDIST, XVAR) 


TOT=0.0 
SUBTOT=1.0 


DO 10 I=1,IDIV/2-1 


SUBTOT=XVAR*SUBTOT 
TOT=TOT+SUBTOT 


210 CONTINUE 
SUETOT= (052.2) / (TOT+1) 
IF(IA.LE.IDIV/2) THEN 
IB=IA 


ELSE 











IB=1A- (IDIV/2) 
ENDIF 
IF (IB.EQ.1) THEN 
DELDIST=SUEBTOT 
ELSE 
DO z0 I=1,IB-1 
SUBTOT=XVAR*SUBTOT 
20 CONTINUE 
DELDIST=SUBIC7 
ENDIF 
RETURN 
END 
ceceecceeceecececeeeecececeeccececeeecececeeececcececceeccececece 
eceecececcececceeceeecceececececeeceeecececeecececceccecececcecec 
3 THIS SUBROUTINE IS FOR LINEAR STRETCHING IN THE WAKE X DIR : 
Eos ceecadeviecedscececactuce cccecececancddeceeuceseeeccecaccetece 
ecececceecececceccececeecececeecececececeecececcececececcececcce 
SUBROUTINE STRCH5 (DIA, IDIV,IA,DELDIST, XVAR) 


TOT=0.0 
SUBTOT=1.0 


DO 10 I=1,IDIV-1 


SUBTOT=XVAR*SUBTOT 
TOT=TOT+SUBICT 


10 CONTINUE 
SUBTOT=DIA/ (TOT+1) 
IF(IA.EQ.1) THEN 
DELDIST=SUBTOT 
ELSE 
DO 20 I#1],IA-1 
SUBTOT=XVAR*SUETOT 


20 CONTINUE 

















DELDIST=SUBTOT 
ENDIF 
RETURN 


END 


/eccecececececceececcececcececececececececcecececccccececccccccce 
cececececcececcecececceccecccecececececceeccececececececececccecccc 
c Cc 
C THIS SUBROUTINE IS FOR LINEAR STRETCHING IN THE Y DIRECTION C 
Cc Cc 
ccececececececececececececcccececceececececeeccecececececcececececce 
ececececccecceeccecececcecccececceceeecececececececececcececececce 

SUBROUTINE STRCH6(DIA,IDIV,IA, DELDIST, XVAR1,XVAR2,INUM, IY) 

TOT=0.0 

SUBTOT=1.0 

XVAR= ((IY-1.0)/ (INUM-1.0) ) *(XVAR2-XVAR1)+XVAR1 

DO 10 I=1,IDIV-1 


SUBTOT=XVAR*SUBTOT 
TOT=TOT+SUBTOT 


10 CONTINUE 
SUBTOT=DIA/ (TOT+2) 
IB=IDIV-IA- i 
IF (IB.EQ.1) THEN 
DELDIST=SUBTOT 

ELSE 
DO 20 I=1,IB-1 

SUBTOT=XVAR*SUBTOT 
20 CONTINUE 

DELDIST=SUBTOT 

ENDIF 

RETURN 

Eno 


ERAERARRARE NEE RHEE KERR ERAERE RR RRREREAREE REET EREEEERRAREREEREEEREREEEK 
RRR RRRARRRERERERE ERE REERERARREE ERE ERE RE REREREEREREREEREAREEEEREEAER ERED 


* THIS IS A SUBROUTINE FOR EXPONENTIAL STRETCHING IN THE Y-DIRECTION * 
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RHEE REE RAKE EEREKEEEKEKERRERERERREEREKKEREKKEEEKKEEKEAKEEEEEEKEK 
ERR REKEKEEREERRREREEEKKEAKEEEREEERREKREKEKEEREKEKEKEKKEKEKEKEKEEEKEK 


SUBROUTINE STRCH? (DIA, IDIV, 1A, DELDIST, GUN, CON2,EX1,EX2, KX, KA) 
CON=(KA-1.0) /(KX-1.0) * (CON2-CON1)+CON1 
EX=(KA-1.0) / (KX-1.0) * (EX2-EX1) +EX1 

DIA1=CON*DIA**EX 

DELY=DIA1/IDIV 


DISTY1=IA*DELY 
DISTY2=(IA-1.0) *DELY 

VAL1= (DISTY1/CON) ** (1.0/EX) 
VAL2= (DISTY2/CON) ** (1.0/EX) 
DELDIST=VAL1-VAL2 


RETURN 
END 


RRR IR RI IT RR TRI TT TR TOR TTR TIRE TOR IT OK tok 
RRR RRR RARER KERR RRR REE REE RRR REE REE EERE RR KER ERE EEER ER IK 
* THIS IS A SUBROUTINE FOR SINEUSOIDALLY STRETCHING THE Y-DIRECTION * 
RRR RIT TR IK TOR KR TIKIT TTR TI TTI TOIT SOI ITT IRR ik 
BIRT RIT I RIT RI RETR RRR EER RE EIR EEK EEK EERE ER ER 


SUBROUTINE STRCH8 (DIA, IDIV, 1A, DELDIST,C1,C2,E1,E2,KX, KA) 


CON= (KA-1.0) / (KX~1.0) *(C2-C1)+C1 
EX=(KA-1.0)/(KX-1.0) *(E2-£1)+E1 
RAD45=.785398163 


DELY=RAD45/IDIV 
DISTY1=IA*DELY+RAD45 
DISTY2=(IA-1.0) *DELY+RAD45 
VALI1=CON*SIN (DISTY1) **EX 
VAL2=CON*SIN (DISTY2) **EX 
DIST=VAL1-VAL2 

DELDIST=DIST*DIA/ (1.0-SIN(RAD45) ) 


RETURN 


END 


RAEKKEREKEEREEKEEKREKREEREREKEREREREKRREKREREAREAEERREERERRAEARERKAEEEEEEEEE 
RHERRKREKEREEEREREERERKEKRREEREEEKEREREKEEREREEEEKEREKEKREREREKEKEREERAKEEKERKEE 
* THIS IS A SUBROUTINE FOR SINEUSOIDALLY STRETCHING THE Y-DIRECTION * 
BRERA REE RRR REREEREEERREREEREREREREEEEREEEEEREAKEEKHAEEKE 
RRR RERRREREE EERE REE EREERERAEREEREKERREEREEKKAKEKEEKEKEE 


SUBROUTINE STRCH9 (DIA, IDIV, IA, DELDIST,C1,C2,E1,E2,KX, KA) 
CON=(KA-1.0)/(K¥=-1.0)*(C2-C1)+¢1 


EX=(KA~-1.0)/(K%-1.0) *(E2-E1)+El 
PI=4. O*ATAN (1.0) 
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RAD45=PI/4.0 
RAD60=PI1/3.0 
RAD90=PI/2.0 


RADSTP=(RAD90-RAD60) / (KX-1) 


DELY= ( ( (RADSTP* (KA-1) ) +RAD60) -RAD45) / (IDIV) 
DISTY1=(IA-1) *DELY+RAD45 

DISTY2=(IA) *DELY+RAD45 

VAL1=CON*SIN (DISTY1) **EX 

VAL2=CON*SIN (DISTY2) **EX 

DIST=VAL2-VAL1 
DELDIST=DIST*DIA/ (SIN (RAD45+DELY*IDIV) -SIN(RAD45) ) 


RETURN 


END 











APPENDIX F. ADDITIONAL SOURCE CODE 
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aaaqaaa 








THIS PROGRAM REMOVES THE FIRST THREE POINTS OF THE HEMISPHERE, BRINGS 
THE APEX TO A POINT, RENUMBERS THE GRID POINTS AND FINALLY DOUBLE THE 
THICKNESS OF THE SURFACE GRID IN THE Z-DIRECTION 


DIMENSION X(140,240,60) ,¥(140,240,60) ,2(140,240,60), 
* XX(140,240,60) , ¥¥(140,240,60) ,22(140, 240,60) 
READ(3) II,JJ,KK 
READ(3) (((X(J,K,L) ,J#1,II) ,K=1,J3J) ,L=1,KK), 
* (((¥(0,K,L) ,J=1,11) ,K=1,37),L=1,KX), 
* (((2(5,K,L) ,J=1,1I) ,K=1,39) ,L=1, KK) 
DO 10 I=1,II-3 
DO 10 J=1,JI 
IF(I.EQ.1) THEN 
XX(1,7,1)=0.0 
YY(1,7,1)=0.0 
Z2Z(1,73,1)=0.0 
ELSE 
XX(I,0,1)=X(I+3,7,1) 
YY¥(I,7,1)=Y¥(I+3,9,1) 
2Z2(I,0,1)=Z(I+3,0,1)*2.0 
ENDIF 
10 CONTINUE 
REWIND 3 


WRITE(3) II-3,37,1 
WRITE(3) ((XX(J,K,1),J=1,II-3) ,K=1,33), 


* ((¥¥(J,K,1) ,J=1,II-3) ,K=1,JJ), 
* ((Z2(3,K,1) ,J=1,II-3) ,Kel,JJ) 
STOP 
END 
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Cc 

C THIS PROGRAM AGAIN DELETES THREE GRID POINTS FROM THE NOSE OF THE 

C SURFACE GRID. THIS PARTICULAR PROGRAM THEN CLOSES THE APEX OF THE 

C SURFACE GRID TO A POINT AND FINALLY RENUMBERS THE GRID FOR FIELD GRID 
C GENERATION 

Cc 


DIMENSION X(150,240,60) ,¥(150, 240,60) ,2(150,240,60), 
* XX (150, 240,60) ,¥Y(150,240,60) ,22(150, 240,60) 
READ(3) II,JJ,KK 
READ(3) (((X(J,K,L) ,J=1,II) ,K=1,JJ) ,L=1,KK), 
* (((¥(9,K,1L) ,J=1,11) ,K=1,33) ,L=1,KK), 
* (((2(0,K,L) ,J=1,11) ,K=1,33) ,L=1, KK) 
bo 10 I=1,II-3 
DO 10 J=1,37 
IF(I.EQ.1) THEN 
XX(1,J7,1)=0.0 
¥Y(1,3,1)=0.0 
Z2(1,7,1)=0.0 
ELSE 
XX(I,J,1)=X(I+3,J5,1) 
¥Y¥(I,d,1)=¥(I+3,5,1) 
22(I,7,1)=Z(1I+3,J,1) 
ENDIF 
10 CONTINUE 
REWIND 3 
WRITE(3) II-3,37,1 
WRITE(3) ((XX(J,K,1) ,J=1,II-3) ,K=1,JJ), 


* ((YY(J3,K,1) ,J=1,TI-3) ,K=1,J39), 
- ((2Z2(3,K,1) ,J=1,II-3) ,Ke1,JJ) 
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VALUE BY DOUBLING THE THICKNESS 


c 
¢e 
Cc 
c 


THIS PROGRAM READS THE FIELD GRID AND JUST CHANGES THE Z-DIRECTION 


DIMENSION X(150,240,60) ,¥(150,240,60) ,2(150,240,60), 
* XX(150,240,60) ,YY(150,240,60) ,2Z2(150,240,60) 


READ(3) II,J3,KK 


READ(3) (((X(J,K,L) ,J=1,II) ,K=1,JJ) ,L=1,XK), 
ig (((¥(9,K,L),3=1,II) ,K=1,3J) ,L=1,kK), 
* (((2(5,K,L) ,J=1,II) ,K=1,3J) ,L=1,KK) 


bdo 10 I=1,II 
DO 10 J=1,d7 
IF(I.EQ.1) THEN 
XX(1,5,1)=X(I,J5,1) 
YY(1,9,1)=0.0 
22(1,9,1)=0.0 
ELSE 
XX(I,J,1)=X(I,J,1) 
YY(I,J,1)=Y¥(1I,J3,1) 
22(1,9,1)=2(I,73,1)*2.0 
ENDIF 
10 CONTINUE 
REWIND 3 
WRITE(3) II,JJ,1 
WRITE(3) ((XX(J,K,1),J=1,II) ,K=1,73), 
* 


((¥¥(J,K,1),J=#1,II),Ke1,JJ), 
* ((22(5,K,1),J=1,11I) ,K=1,JJ) 
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c 

C THIS PROGRAM ALSO DELETES THE FIRST THREE POINTS OF THE HEMISPHERE ON 
C THE SURFACE GRID. BUT THIS PROGRAM ALSO DOUBLES THE THICKNESS OF THE 
C SURFACE GRID IN THE Z-DIRECTION ONLY 

c 


DIMENSION X(150,240,6C) ,¥(150,240,60) ,2(150,240,60), 
* XX(250,240,60) , Y¥(150,240,60) ,22(150,240, 60) 


READ(3) II,JJ,KK 
READ(3) (((X(9,K,1L) ,J=1,II) ,K=1,J3) ,L=1,KK), 
* (((¥(J,K,L) ,J=1,II) ,K=1,J7) ,L=1,KK), 
* (((Z(3,K,L) ,J=1,II) ,K=1,30) , L=1, Kk) 
DO 10 I=1,II-10 
DO 10 J=1,Ic 
IF(I.EQ.1) THEN 
XX(1,J0,1)=X(I+10,7,1) 
YY(1,7,1)=0.0 * 
22(1,0,1)=0.0 
ELSE 
XX(I,J0,1)=X(I+10,7,1) 
YY(I,0,1)=¥(I+10,39,1) 
22(I,0,1)=Z(I+10,7,1)*2.0 
ENDIF 
10 CONTINUE 
REWIND 3 
WRITE(3) II-10,J3,1 
WRITE(3) ((XX(J,K,1) ,J=1,II-10) ,K=1,JJ), 


* ((¥Y¥(0,K,1) ,J=1,II-10) ,K=1,37), 
* ((Z2Z(3,K,1) ,dJ=1,II-10) ,K=1,J7) 
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THIS PROGRAM READS THE FINISHED SURFACE GRID AND DELETES THE FIRST 
THREE POINTS THAT WERE CONSTRUCTED BY THE HEMISPHERE AND THEN 
RENUMBERS THE GRID FOR USE BY THE FIELD GRID GENERATOR 


Cc 
c 
c 
c 
c 


DIMENSION X(120,240,6G6) ,¥(120,240,65} ,Z2(120,240,60), 
* XX(120,240,60) , YY{220,240,60) ,Z2Z(120, 240,60) 


READ(3) II,JI,KK 
READ(3) (((X(J,K,L) ,J=1,II) ,K=1,33) ,L=1,KK), 
* CC(¥(9,K,1L) ,J=1,1T) /K=1,J3J) ,L#=1,KK), 
* (((2(3,K,L) ,J=1,I1) ,K=1,33) ,L=1, Kk) 
DO 10 I=1,II-10 
po 10 J=1,33 
IF (I.EQ.1) THEN 
XX(1,J,1)=X(I+10,J,1) 
YY(1,3,1)=0.0 
2Z(1,J,1)=0.0 
ELSE 
XX (I,J,1)=X(I+10,J,2) 
| ¥¥(1,3,L)=¥(I+1G6,J,1) 
| Z2(1,J,1)=Z(I+10,J,1) 


ENDIF 


10 CONTINUE 
REWIND 3 
WRITE(3) II-10,37,1 
WRITE(3) ((XX(J,K,1) ,J=1,II-10) ,Kel,JJ), 


* ((¥Y(3,K,1) ,J=1,II-10) ,K#1,J7), 
* ((22(3,K,1) ,J=1,II-10) ,Ke1,J37) 
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c 
C THIS IS A PROGRAM THAT MANUALLY ADJUSTS THE FIELD GRID NOSE REGION 
C FOR THE CYLINDRICAL GRID. IT REDISTRIBUTES THE X VALUE ALONG THE 
C SINGULARITY OF THE NOSE. {IT ALSO RESIGNS THE Y AND Z VALUES FOR THE 
C FIELD GRID. 
c 
DIMENSION X(117,240,64) ,¥(117,240,64) ,2(117,240,64), 
* XX(14,240,64) ,YY(14,240,64) ,2Z(14,240,64), 
* XXX(130,240,64) ,YYY(130,240,64) ,2ZZ(130,240, 64) 
OPEN (UNIT=10, FILE='ddwnos. in', STATUS='OLD') 
READ(10,*) DIA,XVAR 
READ(24) II,JJ,KK 
READ(24) (((X(J,K,L) ,J=1,I11) ,K=1,J7) ,L=1,KK), 
* (((¥(J,K,L) ,J=1,11) ,K=1,JJ) ,L=1,KK), 
* (((2(3,K,L) ,J=1,II) ,K=1,J7) ,L=1, KK) 
DO 10 I=1,J7 
DO 10 J=1,kK 
XX(1,I,7)=X(1,I,J) 
Y¥(1,1,7)=Y(1,I,J) 
ZZ(1,1,5)=2(1,1,9) 
10 CONTINUE 
DO 13 I=2,14 
CALL STRCH(DIA,13,I-1,DELDIST,YVAR) 
DO 13 J=1,ag 
DO 13 K=1,KK 
IF(I.EQ.2) THEN 
XX(I,J,K)=DELDIST 


ELSE 


XX(I,0,K)=®XX(I-1,7,K)+DELDIST 


ENDIF 
13 CONTINUE 


DO ll I=2,14 
DO 11 J=1,93 
DO 11 K=1,KK 
IF(K.EQ.1) THEN 


YY¥(I,5,K)=0.0 
22(I,35,K)=#0.0 











ELSE 


YY(1I,J,K)=¥(1,3,K)-Y(1,7,1) 
22(1,0,K)=2(1,3,K)-2(1,9,1) 


ENDIF 
11 CONTINUE 
DO 17 I=1,130 
DO 17 J=1,JI 
DO 17 K=1,KK 
IF(I.LE.13) THEN 
XXX(I,J,K)=XX(15-1,7,K) 
YYY(I,J,K)=¥Y¥(1I,J3,K) 
222Z(1,d,K)=Z2(1,J,K) 
ELSE 
XXX(1I,J,K)=X(I-13,7,K) 
YYY(I,J,K)=¥(I-13,J,K) 
Z22(1,7,K)=2(I-13,3,K) 


ENDIF 
17 CONTINUE 


REWIND 34 


WRITE(34) 130,J3,KK 

WRITE(34) (((XXX(J,K,L) ,J=1,130) ,K=1,JJ),L=1,KK), 
(((¥YY(3,K,L) ,J=1,130) ,K=1,3J) ,L=1,KK), 
(((Z22(5,K,L) ,J=1,130) ,K=1,JJ) , L=1, KK) 


eccececeecceccecceccccceccceccceccecececececeececcceecceccececcecccccececccece 
ececcececececcececcecceccceccceccceceececcecececececceeececcecccccecccc 
C THIS IS THE SUBROUTINE THAT REDISTRIBUTES THE GRID POINTS USING A C 
C LINEAR TYPE DISTRIBUTION. c 
eececceceeceececeeceeecececceceeececeeeececeeeeeeeeeeeeceeeceeeeeeceece 
ecceeecececeececceccceeceececcececeececececeececececececcecccccececcece 


SUBROUTINE STRCH(DIA,IDIV,1IA,DELDIST, XVAR) 


TOT=0.0 
SUBTOT=1.0 


DO 10 I=1,IDIV-1 


SUBTOT=XVAR*SULTOT 
TOT=TOT+SUBTOT 
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CONTINUE 


SUBTOT=DIA/ (TOT+1) 


IF(IA.EQ.1) THEN 
DELDIST=0.0-SUBTOT 
ELSE 
DO 20 I=1,IA-1 
SUBTOT=XVAR*SUBTOT 
CONTINUE 
DELDIST=0.0-SUBTOT 
ENDIF 
RET™RN 


END 














Cc 

C THIS IS A PROGRAM THE READS THE FINISHED SURFACE GRID AND CAN 

C BUILD A FILE OF ANY PARTICULAR YZ-PLANE CROSS SECTION THAT IS DESIRED. 
C THE DESIRED PLANE IS CHOSEN BY CHANGING THE "I" VARIABLE. 

Cc 


DIMENSION X(117,240,60) ,Y(117,240,60) ,2(117,240,60), 
* XX(1,240,60) ,YY(1,240,60) ,22(1,240, 60) 


READ(24) II,JJ,KK 
READ(24) (((X(J,K,L) ,J=1,II),K=1,JJ) ,L=1,KK), 
. (((¥(3,K,L) ,J=1,II),K=1,3) ,L=1,XKK), 
* (((2(3,K,L) ,J=1,1I) ,K=1,3J7) ,L=1,KK) 
I=89 
DO 10 J=1,97 
DO 10 K=1,KK 
XX(1,7,K)=X(I,7,K) 
YY(1,9,K)=¥(1I,J,X) 
22(1,3,K)=2Z(1,9,K) 
10 CONTINUE 
REWIND 37 


WRITE (37) 1,J3,KK 
WRITE(37) ((XX(1,K,L) ,K=1,55) ,L=1,KK), 


* ((¥¥(1,¥,1),K=1,JJ),L=1,KK), 
* ((Z2(1,K,2) ,K=1,J2) ,L=1, KK) 
STOP 
END 
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